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Abstract:  Nonlinear  adaptive  wavelets  are  developed  for  compact  data  representation  and 
efficient  flow  solution  algorithms.  To  flow  solution  algorithms,  proposed  new  adaptive  wavelets 
enhance  the  computational  efficiency  as  well  as  preserve  the  numerical  accuracy  in  steady  and 
unsteady  flow  computations.  Theses  advantageous  features  of  the  wavelet  are  confirmed  by  the 
numerical  tests  of  two-dimensional  airfoil  problem  and  unsteady  shock-vortex  interaction 
problem.  An  adaptive  wavelet  method  with  higher  order  of  accuracy  is  also  proposed  to  allow 
more  accurate  flow  computation.  For  the  high  order  accuracy  of  a  solution,  higher  order  of 
interpolating  polynomial  is  utilized  in  wavelet  decomposition  process.  This  high  order  adaptive 
wavelet  method  was  successfully  applied  to  one-dimensional  shock-sine  interaction  problem, 
two-dimensional  shock-vortex  interaction  problem,  isentropic  vortex  problem,  and  viscous  shock 
tube  problem.  Through  these  applications,  the  computational  efficiency  is  substantially  enhanced 
while  higher  order  numerical  accuracy  of  a  CFD  solver  is  successfully  preserved. 

Introduction:  In  various  aerospace  transportation  and  exploration  systems,  the  prediction 
capability  of  flow  including  supersonic/hypersonic  phenomena  is  highly  important  in  the  design 
of  aerospace  vehicles.  For  the  improvement  of  the  capability,  high  fidelity  simulation  of 
Computational  Fluid  Dynamics  (CFD)  is  required  that  entails  the  development  of  high  order 
method  of  discretized  equations  and  large  number  of  grid  point  clustering.  However,  current 
level  of  technologies  often  falls  short  of  the  need  to  meet  both  accuracy  and  efficiency  at  the 
same  time.  Accuracy  order  increment  and  dense  grid  system,  as  an  effort  to  unveil  complex 
physical  flow  phenomena,  eventually  pay  back  with  substantial  increase  of  computation  time. 
Hence,  it  is  important  to  accurately  compute  the  flow  phenomena  and  at  the  same  time  it  should 
efficiently  alleviate  the  computational  workload. 

For  the  efficient  and  accurate  computations  of  the  flow  phenomena,  there  have  been  many 
attempts  for  the  advancement  of  CFD  algorithms:  one  of  the  representative  examples  is  the 
higher  order  spatial  interpolation  scheme  with  the  limiting  process.  The  main  purpose  of  higher 
order  scheme  is  the  reduction  of  the  number  of  grid  points  without  the  loss  of  the  accuracy  of 
CFD  calculations.  However,  high  order  computations  always  involve  the  unnecessary  numerical 
oscillations  near  a  discontinuity  whenever  shock  waves  are  present.  Because  these  numerical 
oscillations  deteriorate  the  stability  and  convergence  properties  of  the  CFD  computations,  there 
is  an  additional  treatment  of  limiting  algorithms  for  oscillation  removal  [1-3].  However, 
computational  cost  of  high  order  accurate  flux  evaluation  is  generally  expensive,  which 
decreases  the  computational  efficiency. 

The  adaptability  of  CFD  data  can  allow  the  possibility  to  overcome  this  defect;  generally  in 
the  CFD  solution,  the  major  parts  of  computational  domain  contain  smooth  flow  pattern  which 
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can  be  simply  and  accurately  interpolated  by  using  the  values  in  neighboring  cells.  A  dense  grid 
is  only  needed  in  the  rapidly  changing  regions  such  as  shock  waves,  vortex  and  boundary  layer, 
etc.  For  that  reason  several  types  of  adaptive  methods,  especially,  the  adaptive  wavelet  methods 
have  been  studied  as  a  way  of  improving  the  computational  efficiency  [4-6]. 

There  are  two  ways  to  apply  wavelets  in  performing  the  flow  analysis.  According  to  the 
classification  by  Vasilyev  and  Kevlahan  [8],  wavelet-based  numerical  methods  for  solution 
algorithm  can  be  classified  as  adaptive  wavelet  Galerkin  methods  (AWGM)  [9-12]  and  adaptive 
wavelet  collocation  methods  (AWCM)  [4-8].  The  major  difference  between  these  two 
approaches  is  as  follows.  According  to  the  research  by  Bacry,  Mallat  and  Papanicolaou  [10]  or 
Holmstrom  and  Walden  [11],  AWGM  solves  PDE  problems  in  a  wavelet  coefficient  space.  It  is 
mentioned  that  nonlinear  operators,  such  as  multiplication,  are  too  computationally  expensive 
when  done  directly  in  a  wavelet  basis.  To  resolve  this  problem,  Beylkin  and  Keiser  [12]  alternate 
the  computing  space  between  the  physical  and  the  wavelet  domain  in  order  to  compute  the 
operators.  In  other  words,  multiplication  is  performed  in  physical  space  and  differentiation  is 
performed  in  wavelet  space. 

On  the  other  hand,  dealing  with  nonlinearity  and  general  boundary  conditions  is  not 
complicated  in  AWCM  because  all  computations  are  performed  on  the  physical  space.  One  can 
directly  link  every  wavelet  coefficients  with  the  collocation  points  on  the  physical  domain.  In 
AWCM  studies,  Harten  presents  an  adaptive  multi-resolution  scheme  for  computing  the 
discontinuous  solutions  of  hyperbolic  PDEs.  [4]  Holmstrom  makes  the  algorithm  for  organizing 
an  adaptive  dataset  with  an  interpolating  wavelet  transformation  [5].  Sjogreen  uses  a 
multi-resolution  scheme  to  solve  compressible  Euler  equations  and  also  utilizes  an  interpolating 
wavelet  transformation  to  construct  the  adaptive  dataset  [6].  He  concludes  that  the 
multi-resolution  method  can  yield  considerable  gains  in  efficiency  by  using  computationally 
expensive  method  with  a  large  number  of  grids. 

However,  the  accuracy  and  efficiency  of  an  interpolating  wavelet  transformation  are  very 
sensitive  to  threshold  values  during  multi-resolution  process.  In  this  research,  numerical  errors  of 
which  the  extent  reaches  a  magnitude  equal  to  that  of  the  threshold  value  inevitably  occur  at  the 
excluded  cells  of  an  adaptive  dataset  and  they  eventually  prevent  the  solutions  from  converging. 
They  are  transferred  to  the  entire  computational  domain  through  a  time  integration  process, 
disturb  the  stable  organization  of  an  adaptive  dataset,  decrease  the  compression  ratio  and  finally 
deteriorate  the  computational  efficiency. 

The  objective  of  the  study  is  the  improvement  of  computational  efficiency  of  CFD 
computation  with  preserving  the  spatial  and  temporal  accuracies  of  conventional  CFD  solvers.  In 
order  to  meet  these  objectives,  improved  adaptive  wavelet  method  includes  the  modified 
threshold  method,  residual  interpolation,  and  restriction  technique  as  follows. 

First,  the  modified  threshold  method  consists  of  the  modified  threshold  value  and  the 
stabilization  process.  The  modified  threshold  value  considers  grid  spacing  and  time  step  for 
retaining  the  original  spatial  and  temporal  accuracy  of  conventional  solvers  by  automatically 
adjusting  the  threshold  value.  And  stabilization  process  is  involved  in  stable  construction  of  an 
adaptive  dataset  by  obstructing  the  additional  0{s)  error  due  to  threshold  is  transferred  to  the 
adjacent  computational  domain  during  the  time  integration  process.  Therefore,  it  can  contribute 
the  enhancement  of  the  compression  ratio  of  the  current  method. 

Secondly,  general  adaptive  wavelet  procedure  is  modified  by  adopting  residual  interpolation  at 
the  n  time  step,  not  the  n  + 1  time  step.  A  residual  based  interpolation  concept  was  introduced 
by  Chiavassa  and  Donat  for  point  value  setting  [13].  In  this  paper,  this  method  is  modified  for 
cell  value  setting  and  time  integration  is  performed  on  the  entire  domain,  i.e.,  residual  values  at 
excluded  cells  in  an  adaptive  dataset  are  interpolated  by  using  the  calculated  residual  values  at 
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included  cells  in  an  adaptive  dataset.  Because  the  interpolation  is  performed  at  the  n  time  step, 
there  is  no  need  to  add  several  adjacent  cells  to  an  adaptive  dataset.  In  addition,  this  approach 
can  be  applied  to  implicit  time  integration  methods  because  all  residual  values  are  computed. 

Thirdly,  restriction  technique  is  applied  after  the  time  integration  on  the  entire  domain.  In  this 
method,  if  the  flow  variations  in  time  are  small  enough  at  the  excluded  cells  in  an  adaptive 
dataset  when  compared  with  the  order  of  the  threshold  value,  these  variations  are  discarded  and 
restricted  by  multiplying  the  weighting  factor.  This  process  can  especially  contribute  the 
convergence  acceleration  of  steady  state  flow  simulations. 

Throughout  these  processes,  spatial  and  temporal  accuracies  of  a  conventional  solver  are 
preserved  with  a  substantial  reduction  in  computing  time.  Various  steady  and  unsteady  flow 
simulations  are  conducted  to  verify  the  enhancement  of  both  efficiency  and  accuracy  of  the 
proposed  method. 

This  report  is  organized  in  the  following  order.  After  the  introduction,  development  of 
nonlinear  wavelet  is  described  as  a  preliminary  background  in  section  1  of  Experiment.  In 
section  2  and  3  of  Experiment,  we  describe  the  implementation  of  the  new  adaptive  wavelet 
method  for  a  reference  CFD  solver  and  development  of  adaptive  wavelet  method  based  high 
order  scheme.  In  section  1,  2,  and  3  of  Results  and  Discussion,  several  numerical  simulations 
using  adaptive  wavelet  methods  are  performed  and  a  conclusion  is  drawn. 


Experiment: 

1.  Development  of  nonlinear  wavelets 

The  phenomena  of  general  interest  in  CFD,  consist  of  several  flow  features  such  as  Shock, 
Vortex,  Boundary  layer,  Expansion  fan,  etc.  and  those  flow  features  are  depicted  as  follows. 


1  A  ExptJUttUI  WdW 


In  the  analysis  and  representation  of  these  data,  the  wavelets  have  attractive  features  as  locality 
and  compactness.  Through  these  features,  the  wavelets  perform  properly  the  role  of  filtering  that 
classifies  primary  into  a  region  of  rapidly  changing  such  as  shock,  vortex,  boundary  layers,  etc. 
and  a  smooth  region  of  mild  changing  in  flow  properties. 


Let  us  assume  a  dyadic  grid  set  as  shown  in  Eq.  (1-1). 


Vl  =  {%i,k  G  ^  :  xi,k  1  k,  k  e  Z},  /  e  Z  . 


(1-1) 


The  key  idea  of  the  adaptive  wavelet  method  is  that  smoothness  of  flow  pattern  can  be  easily 
determined  by  the  magnitude  of  difference  value  between  original  function  value  and 
approximated  value  based  on  neighboring  cells.  That  is,  the  original  value  can  be  accurately 
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approximated  by  the  values  at  neighboring  cells  in  the  smooth  region  but  not  vice  versa  in 
rapidly  changing  region.  Here,  interpolating  polynomial  is  used  at  the  odd  numbered  cells  for  the 
approximation  of  original  function  values  and  then  the  dyadic  dataset  is  decomposed  as  Eq. 
(1-2). 


fl+\,2k  ~  f l, k 
fl+l,2k+l  =  ^l+l,2k+l(f. 


l,k-—+\  ’ 
2 


'  ’  fl,k>  fl, 


k+\> 


2 


(1-2) 


where  P ^  is  the  order  of  interpolating  polynomial;  the  4th  order  of  interpolating  polynomial  is 
shown  in  Eq.  (1-3). 


19  9  1 

f l+\,2k+\  =  ~~^fl.k-\  +Yk^l’k  +  ~l6^!’k+'  ~  ~l6^l’k+1  ' 


After  the  approximation,  the  difference  can  be  computed  as  Eq.  (1-4). 

d )tk  =  fl+i,2k+i  ~fl+\,2k+i’  VA:  e  Z  . 


(1-3) 


(1-4) 


In  smooth  region,  the  order  of  difference  value  is  very  small  and  the  corresponding  cell  can  be 
partially  rejected  in  the  whole  dataset;  an  adaptive  dataset  is  acquired  according  to  the  local 
features  of  solution.  And  the  maximum  error  is  bounded  within  the  order  of  threshold  value  as 
Eq.  (1-5). 


<0{£). 


(1-5) 


2.  Flow  solution  algorithm  using  adaptive  wavelets  for  CFD  solution  algorithm 

The  two-dimensional  Euler  equations  are  used  as  the  governing  equations  of  the  flow 
problems.  The  two-dimensional  Euler  equations  are  written  as  Eq.  (2-1). 


dQ  dE  dF 

—  H - H - 

dt  8x  dy 


(2-1) 


P 

pu 

pv 

where  Q  = 

pu 

pu2  +  p 

puv 

,  E  = 

,  F  = 

pv 

puv 

pv  +  p 

j Pet_ 

(pp,  +  p)u_ 

(PP,  +  P)v_ 

Z  1 

-t  a  1  ,  2  2x 

and  ct  — - 1 — (u  +  v  ) . 

r(r- 1)  2V 

All  the  properties  and  governing  equations  are  non-dimensionalized. 

By  the  generalized  coordinate  transformation,  Eq.  (2-1)  is  rewritten  as  Eq.  (2-2). 
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where 


Q  =^j,  E  =  -y[£,  +  E,XE  ,+ZyF] ,  F  =^j[rjt  +  r/xE  ,+rjyF]. 

To  improve  the  computational  efficiency  of  steady/unsteady  flow  problems,  the  general 
adaptive  wavelet  procedure  is  changed  by  adopting  residual  interpolation  at  the  n  time  step  and 
time  integration  is  performed  on  the  entire  domain.  The  procedure  of  selecting  and  adding 
additional  cells  to  an  adaptive  dataset  according  to  the  wavelet  resolution  levels  and  positions  is 
excluded.  And  modified  threshold  method  is  applied  in  order  to  maintain  the  spatial  and 
temporal  accuracies  of  reference  schemes.  Also,  a  stabilization  method  is  applied  again  to 
unsteady  problems  to  stably  construct  an  adaptive  dataset  while  the  compressibility  of  the 
wavelet  method  is  enhanced.  Especially  in  steady  state  problems,  the  restriction  technique  is 
applied  for  the  acceleration  of  the  convergence.  The  flowchart  in  Fig.  2-1  shows  the  overall 
implementation  of  the  new  adaptive  wavelet  method  in  a  reference  solver.  This  implementation 
is  comprised  of  the  following  sequential  steps. 


■  Modified  Adaptive  Wavelet  Method  - 


Fig.  2-1  Overall  procedure  of  flow  simulation  with  the  new  adaptive  wavelet  method 

2-1.  Construction  of  an  adaptive  dataset 

The  set  of  dyadic  grids  in  two-dimensions  can  be  presented  as  shown  in  Fig.  2-2  Positions  of 
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the  symbol  O  are  even  numbered  grids  and  the  other  positions  of  symbols  □,  A,  and  X  are 
odd  numbered  cells.  The  4th  order  interpolating  polynomials  are  shown  in  Eq.  (2-3). 
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Fig.  2-2  Example  of  a  two-dimensional  dyadic  grid  set. 


□  : 
A: 

X  : 


eUj  ~a. 

av, = ~q:,-2+^q-,+^q'j„  , 

a+ij+l  =0.5x(- —Q^2j_2  +  —  Q"j  +  —  Q! '+2J+2  |6«;j:.j) 


(2-3) 


Especially  in  unsteady  problems,  a  tiny  change  of  flow  may  grow  into  a  large  change  with 
time.  Thus,  even  a  tiny  variation  has  to  be  considered.  In  this  case,  we  used  the  6th  order  of 
interpolating  polynomial  and  present  the  two  dimensional  formulations  as  shown  in  Eq.  (2-4). 


□  :  ar..j  =^2"u -|A  A®*  ~§< 

a:  p-4) 
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X  : 


aVi  =  x  VQiJ- 4  +  32r+2,y-4  +  20- 2>y_2  -  270”,._2  -  21Q"+2  j_2  +  20/U,,_2 

+  ^Qi-4,j  ~  VQUj  +  174 0",  + 1 740^.2 >7-  -  270”  4,,  +  30"  6>y 
+  30- 4J+2  -  270”_2y-+2  + 1740”,+2  + 1740- 2>>+2  -  27Q"+4j+2  +  30”  6J+2 
+  20”  2J+4  -  27 Q"j+4  -  27Q"+2j+4  +  2Q"+4j+4  +  3 0”,.+6  +  30/(2J+6) 


Here,  the  formulation  at  the  X  position  in  Eq.  (2-3),  (2-4)  uses  24  cells  for  the  6th  order 
interpolating  polynomial  and  is  very  complicated.  Instead  of  this  equation,  we  formulated  Eq. 
(2-5)  using  12  cells  as  follows. 


Qmj+i  =  x  (30"-4 j-4  -  250”_2 j_2  + 1 5O0”y-  + 1  50Q"+2j+2 

-  250”  4J+4  +  30”  6J+6)  +  A  x  (3Q?_4J+6  -  25QUj+4 
+ 1500”, 2  + 1500”  ,  - 250”+4  ,  2  +  30"  ,  4) 


(2-5) 


0 


8  10  12 


Level  3 


X 


Level  2 

Level  1 

o 

□ 

z 

O 

X 

O 

X 

O 

o| 

O 

X 

o 

X 

O 

X 

o 

X 

Ol 

X 

o 

i  =  2^  ^  Cells  :  Always  included  in  a  dataset 


/  =  2l  +  2(/  ^  Cells  :  4th  order  of  interpolation 


Fig.  2-3  Boundary  treatment  in  wavelet  decomposition 

Based  on  many  numerical  simulations,  Eq.  (4-5)  has  much  better  efficiency  than  the  accurate 
formulation  with  little  lose  of  accuracy.  Note  that  some  indices  of  Qij  are  located  outside  of  the 
computation  domain  near  its  boundaries.  In  order  to  avoid  this  situation,  cells  /  =  2(/-1) 
(Loarsest  Level  -  ^  <  ^Finest  Level )  are  always  included  in  the  adaptive  dataset.  At  the  next  odd  numbered 
cells  i  =  21  +  2(/_1) ,  the  4th  order  interpolating  polynomial  is  applied.  The  overall  boundary 
treatment  is  presented  in  Fig.  2-3,  and  the  same  procedure  can  be  applied  to  the  other  side  of 
boundary,  or  the  boundaries  along  the  j -direction. 

Then,  the  wavelet  coefficients  of  each  □ ,  A,  and  X  position  are  found  as  follows: 


□  :  d^lJ=Q^J-Q^ 
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(2-6) 


A: 
X  : 


d"i+l=Qi+l~Q" 


ij+ i  ■ 


d“ 


i+lj+l 


=  GT+u+i-S" 


■i+lj+l 


After  the  wavelet  decomposition,  the  thresholding  process  is  performed.  The  threshold  method 
consists  of  the  modification  of  a  threshold  value  and  the  stabilization  of  an  adaptive  dataset.  By 
this  method,  it  is  possible  to  maintain  the  spatial  and  temporal  accuracies  of  reference  method. 
Also,  due  to  the  stable  construction  of  an  adaptive  dataset,  it  can  enhance  the  compression  ratio 
of  wavelet  method  and  improve  computational  efficiency. 

In  order  to  maintain  the  Ith  order  of  spatial  accuracy  and  mth  order  of  temporal  accuracy  of 
reference  method,  the  modified  threshold  value  is  written  as  Eq.  (2-7). 

s'  =  min[£-,  max(Axz ,  CFLm  ■  Axm )] ,  (2-7) 


where  Ax  is  determined  at  the  finest  resolution  level. 

Based  on  this  modified  threshold  value,  s' ,  We  can  control  the  flag  values  of  the  cells  and 

threshold  a  dataset.  If  d"j  is  larger  than  s' ,  the  flag  value  of  cell  (i,j)  is  determined  as  1  and 
the  cell  is  included  in  a  I(s‘ ')"  dataset;  if  not,  the  flag  value  is  set  as  0  and  the  cell  is  excluded 
from  the  I(s')n  dataset.  Then,  the  dataset  is  adapted  to  the  flow  features  while  maintaining  the 
numerical  accuracy  of  reference  schemes. 

Also,  we  prevent  the  flag  value’s  switching  between  0  and  1  due  to  additional  errors  to 
originate  from  the  organization  of  an  adaptive  dataset  and  realize  the  stable  construction  of  an 
adaptive  dataset  that  leads  to  the  enhancement  of  the  compression  ratio. 

2-2.  Modified  thresholding  method 

Because  |  ^  is  bounded  within  0(s),  °{s)  errors  are  added  to  the  solutions  by  the  wavelet 
transformation.  Then,  each  primitive  variable  or  flux  value  has  0(s)  errors  as  Eq.  (2-8). 


P?j  =  Ph+Oi.e),  u"j  =  u"j  +  O(s)  , 
vF=vlj+0{£),  p"j  =  plj  +  O(s)  , 
E",2=K,2+0(s),  F?I2=F»2+0{£)  . 


(2-8) 


Also,  if  the  flux  values  are  discretized  with  Ith  order  of  spatial  accuracy,  each  primitive 


variable  or  flux  value  has  Ith  order  of  truncation  error  term  as  Eq.  (2-9). 


Pi,j  ~  Pi,),real  +  0(hx'  )  ?  u"j  -  u"j  real  +  0(Ax‘  )  ; 
vi,j  =  V",j,real  +  )  ,  P"j  =  Pineal  +  <2(Ax' )  , 

Kl  =  EWreal  +  <K Ax' )  ,  F”/2  =  F"/2  ,  +  0(Ax' )  . 


(2-9) 


From  the  combination  of  Eqs.  (2-8)  and  (2-9),  the  fluxes  come  to  have  the  following  order  of 
errors  due  to  spatial  discretization  and  thresholding,  as  shown  in  Eq.  (2-10). 

K2  =  El,  2, real  +  0(Ax‘ )  +  O(s),  F?,2  =  F"/2  real  +  0(Ax‘ )  +  0(s)  . 
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(2-10) 


By  mth  order  of  time  integration,  all  terms  of  numerical  errors  are  derived  as  Eq.  (2-1 1). 


At 


2,n+l  _  nn  j?n 

Li  ~  &ij  \F/2,real  ^-1 


Ax 


1/  2, real 


1/  2, real 


F-M2  ,real) 


(2-11) 


+  0( Ax1  ■  At)  +  0(Atm+l)  +  0(s  •  AO 

Therefore,  0(e)  have  to  be  smaller  than  0(Axl)  and  0(Atm)  in  order  not  to  damage  the  Ith 
order  of  spatial  accuracy  and  mth  order  of  temporal  accuracy  of  the  conventional  method.  For  the 
satisfaction  of  this  purpose,  modified  threshold  value  is  developed  as  Eq.  (2-12)  [14,  15]. 


s'  =  min[£,  max(Ax/ ,  CFLm  •  Axm )]  .  (2- 1 2) 

Because  our  target  is  3rd  order  accuracy,  3rd  order  of  spatial  and  temporal  accuracy  of  the 
numerical  schemes  are  used.  Then,  the  modified  threshold  value  in  this  research  is  proposed  as 
Eq.  (2-13).  where  CFL  means  the  constant  of  CFL  conditions. 


s'  =  min^maxi Ax3, CFL3  •  Ax3)]  . 


(2-13) 


Based  on  this  modified  threshold  value,  the  odd  numbered  cells  are  included  in  the  dataset  if 
dij  is  larger  than  s' ;  if  not,  the  cells  are  excluded  from  the  dataset.  Then,  the  dataset  is 
constructed  to  follow  the  local  features  of  the  solution  with  the  conservation  of  the  numerical 
accuracy  of  the  conventional  CFD  schemes. 

2-3.  Flux  evaluation  at  included  cells  in  an  adaptive  dataset 


After  constructing  the  I(e')n  dataset,  the  flag  values  of  the  dataset  determine  whether  the  flux 
values  are  to  be  evaluated  or  not.  If  the  flag  value  is  1,  the  flux  value  is  calculated  by  reference 
schemes  such  as  AUSMPW+  or  Roe’s  FDS  method.  In  this  research,  we  use  the  Finite  Volume 
Method  (FVM);  i.e.  the  numerical  fluxes  E,  F  are  defined  at  these  interfaces.  Thus,  if  the  flag 
value  of  cell  (Uj)  is  1,  flux  evaluations  are  performed  on  all  sides  of  the  cell,  and  residual 
values  at  cell  (hj)  are  calculated.  If  the  flag  value  is  0,  flux  evaluation  is  not  performed,  so  the 
total  computation  time  is  decreased. 

2-4.  Stabilization  of  an  adaptive  dataset 


If  the  solution  changes  after  the  time  integration  step,  the  adaptive  dataset,  I(£')n  also 
changes  in  order  to  be  adapted  to  the  latest  solution.  Independent  of  a  solver,  s’  order  errors  are 
suddenly  imposed  due  to  changes  of  the  adaptive  dataset  and  are  transferred  to  the  surroundings 
throughout  the  time  integration  step.  Thus,  both  the  surroundings  and  wavelet  coefficients  at  the 
surroundings  have  the  s’  order  errors.  The  s'  order  errors  in  the  wavelet  coefficients  obstruct 
the  construction  of  an  adaptive  dataset  adapted  to  the  solution.  Consequently,  the  compression 
ratio  becomes  reduced  and  the  computational  efficiency  of  the  overall  adaptive  wavelet  method 
decreases.  For  the  efficient  construction  of  an  adaptive  dataset,  the  alteration  of  flag  values  due 
to  order  errors  should  be  restricted. 

We  use  algorithm  1  (defined  below)  for  keeping  I(e')n  from  being  switched  regardless  of  the 
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solution’s  properties.  First,  I  compare  I(s’)n  with  I{sf)n  1 ,  in  the  case  where  the  point  was 
already  included  in  I(s')n~l ,  the  point  can  be  excluded  from  I(s')n  only  if  the  wavelet 
coefficient  at  that  point  is  smaller  than  as '  ( oc  <  1)  ?  not  smaller  than  s' .  On  the  other  hand,  in 
the  case  where  the  point  was  already  excluded  from  I(s')n _1 ,  then  the  point  can  be  included  in 
I(s')n  if  the  wavelet  coefficient  is  larger  than  fie*  ( P  >  1 )  but  not  sf . 


2-5.  Switching  of  residual  calculation  method 

Generally,  flow  properties  are  interpolated  after  time  integration  in  previous  wavelet  methods. 
That  is,  flow  properties  are  reconstructed  at  the  n  + 1  time  step  based  on  the  I(s')n .  This 
discrepancy  of  time  step  causes  the  following  numerical  problem:  if  the  flow  changes  rapidly, 
e.g.,  the  shock  or  vortex  propagation,  etc.,  the  neighboring  cells  of  the  cell  that  is  included  in 
I(s')n  become  important  as  well  as  the  cell  in  I(s')n  itself.  If  they  are  not  included  in  lie  y, 
the  interpolation  at  these  neighboring  cells  cannot  guarantee  the  solution  accuracy  and  may 
deteriorate  it  at  the  n  + 1  time  step. 

For  resolving  this  problem  residual  interpolation  process  is  introduced  in  this  study.  Residual 
values  at  excluded  cells  in  an  adaptive  dataset  are  calculated  by  interpolation  polynomials  which 
uses  the  values  at  even  numbered  positions  (O  positions)  in  Fig.  2-2.  The  4th  order  equations  for 
residual  interpolation  are  as  follows: 


□  : 

A: 


R 


■*'+1,7 


Kj* 


19  9  1 

= - R'\  ,  +—R";  +—RL,  , - R, 


16 


li-2J 


16 


16 


li+2J 


16 


J+4J  » 


1  9  9 

-R"_,+  —  R''+  —  R" 


16 


2 


16 


16 


i,j+ 2 


16 


-R" 


i,j+ 4  ’ 


X: 


0-5X^R'-W  +  T6R",,2  +  ^K2,J--^K,J-1) 


(2-14) 


Especially  in  unsteady  problems,  for  representing  the  tiny  variations  of  flow  properties,  the  6th 
order  interpolating  polynomial  is  used  as  shown  in  Eq.  (15). 


□  : 

A: 

X: 


j?n  =  ^  pn  25  n  75 

iJ+l  256  iJ~4  256  iJ~2  128 J 


1ZO 


1 


rx  (3^4,7-2-25^ 


■2,2-2  ' 


DH  _  1 

M’J+l  ~  512 

A  x  (3 R?_4J+6  -  25R"_2j+4  +  150/?,"/+2 


150^,+ 150^2,  j+2 


+  150  R"+2J 


-25R''+4j+4  +  3R’'+6j+6) 
+3  Ku-t) 


-25^  j- 


After  the  construction  of  residual  values  in  the  whole  computational  domain,  time  integration  is 
performed  by  using  conventional  schemes. 
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3.  Development  of  adaptive  wavelet  based  high  order  scheme 


3-1.  High  order  spatial  interpolation  schemes  for  CFD 

High-accuracy  numerical  methods  are  required  to  obtain  converged  solution  as  well  as  to 
capture  correct  flow  physics.  High-fidelity  simulation  calls  for  high-accuracy  algorithm  and/or 
dense  grid  system,  leading  to  substantial  computation  workload.  Thereby,  it  is  very  important  to 
develop  an  efficient  numerical  scheme  that  achieves  high  order  of  accuracy  at  the  same  time 
captures  correct  unsteady  flow  phenomena.  There  are  various  high  order  numerical  methods 
developed  to  date,  such  as  TVD  (Total  Variation  Diminishing)  scheme  [16-18],  ENO/WENO 
(Essentially  Non-Oscillatory  /  Weighted  Essentially  Non-Oscillatory)  scheme  [19,  20],  MLP 
(Multi-dimensional  Limiting  Process)  scheme  [2],  e-MLP  (enhanced  Multi-dimensional  Limiting 
Process)  scheme  [3],  Galerkin  method.  Spectral  and  Spectral  volume  methods. 


Table  3-1.  Hig: 

her  order  schemes  for  CFD 

TVD 

ENO 

MLP 

e-MLP 

Discontinuity  in  one  dimension 

O 

O 

O 

O 

Discontinuity  in  multi-dimension 

X 

X 

O 

O 

Local  extrema 

X 

o 

X 

O 

Those  offer  some  advantageous  features  but  none  of  them  is  universally  better  than  the  others. 
In  TVD,  incorrect  results  such  as  clipping  phenomena  appears  at  local  extrema,  and  generally  not 
easy  to  extend  to  multi-dimensions  since  it  was  developed  in  ID  based  equations.  The  usage  of 
ENO  scheme  can  avoid  clipping  phenomena,  but  undershoot/overshoot  may  appear  in 
multi-dimensional  analysis.  In  addition,  the  use  of  unstructured  mesh  gives  some  difficulty  in  the 
construction  of  an  appropriate  stencil,  leading  to  excessive  memory  load.  MLP  shows  a 
monotonic  and  stable  calculation  of  the  shock  wave  in  multi-dimensional  problems.  However, 
MLP  damps  out  the  flow  and  is  unsuitable  for  local  extrema  in  a  continuous  region  as  with  other 
limiting  functions,  which  decreases  the  accuracy  of  the  solutions.  e-MLP  has  distinguishing  step. 
Through  this  step,  the  computational  domain  is  divided  into  continuous,  linear  discontinuous  and 
nonlinear  discontinuous  regions.  Appropriate  limiting  criteria  are  then  applied  to  each  region. 
This  increases  the  accuracy  of  the  solution  and  makes  it  robust  against  shock  instability.  For  the 
development  of  adaptive  wavelet  method  based  high  order  scheme,  we  adopt  TVD  and  e-MLP 
scheme  respectively. 

3-2.  High  order  accurate  wavelet  decomposition 

If  an  initial  two-dimensional  dyadic  grid  set  with  level  l  is  use,  as  shown  in  Fig.  5-1,  then  O 
cells  are  even  numbered  and  the  others  are  odd  numbered  cells. 
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i-3  i-2  i-1  i  i+1  i+2  i+3  2i-6  2i-4  2i-2  2i  21+2  21+4  2i+6 

j+3 

j+2 

j+1 

j 

j-1 
j-2 
j-3 

Fig.  3-1  Example  of  a  two-dimensional  dyadic  grid  set. 
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(a)  Level  l 


2j+6  oooo ooo ooo oo o 
AXAXAXAXAXAXA 
2j+4  oo  oo  oo  o  ooo  oo  o 
AxAxAxAxAxAxA 

2j+2  ooooooooooooo 

AxAxAxAxAxAxA 

2j  ooooooooooooo 

AxAxAxAxAxAxA 

2j-2  ooooooooooooo 

AxAxAxAxAxAxA 

2j-4  ooooooooooooo 

AxAxAxAxAxAxA 

2j-6  loloblolololololololooo 

(b)  Level  /+ 1 


The  values  at  even  numbered  cells  are  saved  to  cells  in  the  coarser  level  grid.  In  the  other  odd 
numbered  cells,  an  interpolating  polynomial  is  applied  to  approximate  the  values  using  the 
values  at  the  even  numbered  cells,  as  in  Eq.  (3-1).  Where  is  the  interpolating  polynomial, 
and  p  is  an  even  number. 

Wavelet  coefficients  at  odd  numbered  cells  are  then  calculated  as  the  difference  between  the 
solution  values  and  the  approximated  values.  In  order  to  maintain  the  original  numerical 
accuracy  of  the  conventional  CFD  method,  higher  order  accurate  wavelet  transformation  is 
needed.  For  the  selection  of  the  order  of  interpolation  in  the  wavelet  decomposition,  let  us 
consider  the  spatial  discretization  method  with  3rd  order  accuracy.  The  residual  values  in  this 
case  have  the  4th  order  of  the  spatial  error  term  as  given  by  Eq.  (3-2). 


*  Av(e;i  -  E’_ t  +  ot&x1 ))+  A»(f"4  -  F’ +  0(A>.3)j 
=  R":  +  0(Ax2Ay)  +  0(Ay2Ax) 


:R",+0(Ax4) 


12 


From  Eq.  (3-2),  the  5th  order  of  the  interpolating  polynomial  is  sufficiently  accurate  and  is 
recommended  for  the  accurate  reconstruction  of  residuals  in  the  excluded  cells  in  the  adaptive 
dataset.  Especially  in  unsteady  flow  computations,  an  at  least  6th  order  interpolating  polynomial 
is  needed  for  more  precise  calculation  of  the  local  flow  features  in  an  unsteady  problem  [8],  For 
consistent  accuracy,  the  same  order  of  interpolating  polynomials  is  necessary  not  only  for  the 
residual  interpolation,  but  also  for  interpolation  in  wavelet  decomposition.  Therefore,  the  6th 
order  interpolating  polynomial  is  used  for  the  wavelet  decomposition  process  as  well  as  for  the 
residual  interpolation  process.  Eq.  (3-3)  shows  the  details  of  the  sixth-order  wavelet 
decomposition  process. 


O : 
A: 


Qmj  =  ^(30+4,2  -250/l2, j  +  1500”,.  +  15O0”+2J -25Q"+4j  +  3 Q"+6J) 

Z  JO 

Qmj  =  t~t(30/j-4  - 25 QIj_2  +  150 Qlj  +  150 Q?J+2  -  25 0”y+4  +  3 Q"j+6) 

ZJO 

0, ”+1,2+1  =  X  (300-4  +  30,4-2,7-4  +  20”  2,7-2  -  270^-2  ~  27 0”  2,7-2  +  20”  4,7-2 

+  30O4,7  - 27002,7  + 17400  + 1742,” 2,7  - 270”4J  +  30”  6J 
+  30”  4,7+2  -  27002,7+2  +17400+2  +174002,7+2  -  270O4,7+2  +  30O6,7+2 
+  2002,7+4  -  2700+4  -  27002,7+4  +  20”  4,7+4  +  300+6  +  30”  2,7+6) 


(3-3) 


After  interpolation  of  the  solution  data  at  the  odd  numbered  cells,  the  difference  values,  d"j  of 
each  O  ,  A  and  X  position  are  calculated  as  in  Eq.  (2-4). 


O: 
A: 
x  : 


dUj=QUj~Qi\j 

ao+i-ao+i 

=  0,01,2+1-0,+ 


rl" 

ai,j+ 1 


d 


,+  1,2+1 


,+  1.2+1 


(3-4) 


The  decomposition  process  is  performed  at  the  next  coarseness  level  of  grid  set  and  all  d"j 
are  calculated  with  multi-resolution. 


3-3.  Preservation  of  high  order  accuracy 

Because  \d"j  \x  is  bounded  within  0(e) ,  0(e)  errors  are  added  to  the  solutions  by  the 

wavelet  transformation.  Thus,  each  primitive  variable  or  flux  value  has  0(e)  errors  as  in  Eq. 
(3-5). 


P^Pij+Oie),  u^j=KJ+0(e), 

=  vUj  +  0(e) ,  pIj=pIj+0(£), 

A"  2  =  K  2  +  0(s) ,  F”2  =  f;)2  +  0(8) . 


(3-5) 


Also,  if  the  flux  values  are  discretized  with  the  Ith  order  of  spatial  accuracy,  each  primitive 
variable  or  flux  value  has  a  Ith  order  truncation  error  term,  as  in  Eq.  (3-6). 
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(3-6) 


PiJ  =  Pi,j,real  +  0( Ax')  ,  u?j  =  <;,rea/  +  0( Ax')  , 
v"|  =  VUj,real  +  0(AxZ  )  #  p"j  =  p"jtreat  +  0(Ax‘  )  ? 

^f/2  =  +  0(Ax' )  ,  F,"2  =  Fyllreal  +  0(Axl ) 

Combining  Eq.  (3-5)  and  (3-6),  the  fluxes  have  the  following  order  of  errors  due  to  spatial 
discretization  and  thresholding,  as  shown  in  Eq.  (3-7). 

Kl  =  K2,real  +  <K V  )  +  0(f),  +  0(Ax‘ )  +  0(f)  (3-7) 

Then,  by  the  mth  order  of  time  integration,  all  terms  of  the  numerical  errors  are  derived  as  Eq. 
(3-8). 


Qu  Qu  ~  7 


!2,real  1/2, real  +  ^1/2, real  E_|  /  2, real ) 


’  Ax 

+  0( Ax'  •  Ar)  +  0(Afm+1)  +  0(s  ■  At) 


(3-8) 


Therefore,  0(f)  has  to  be  smaller  than  O(Ax')  and  0(Atm)  in  order  to  avoid  damaging  the 
Ith  order  of  spatial  accuracy  and  the  mth  order  of  temporal  accuracy  of  the  conventional  method. 
To  this  end,  the  modified  threshold  value  is  developed  as  in  Eq.  (3-9)  by  assuming  a  linear 
advection  equation  and  introducing  the  CFL  condition. 

f'  =  min[f,max(Ax/,CFAm -Ax'")]  (3-9) 


Because  our  target  is  third-order  accuracy,  the  third-order  spatial  and  temporal  accuracies  of  the 
numerical  schemes  are  used.  The  modified  threshold  value  in  this  research  is  then  proposed  as 
Eq.  (3-10). 


s’  =  min[f,max( Ax3, CFl}  Ax3)]  (3-10) 

Based  on  this  modified  threshold  value,  s’ ,  the  odd  numbered  cells  are  included  in  the  dataset 
only  if  d"j  is  greater  than  s' .  The  dataset  is  then  constructed  to  follow  the  local  features  of  the 
solution  while  preserving  the  numerical  accuracy  of  the  conventional  CFD  schemes. 

3-4.  Flux  evaluation 

In  view  of  Godunov-type  approach,  the  steps  to  construct  a  numerical  flux  at  a  cell-interface 
usually  consist  of  interpolation  stage  and  evaluation  stage.  It  is  known  that  interpolation  stage  is 
generally  independent  of  evaluation  stage  where  the  local  Riemann  problem  is  solved  at  a  cell 
interface.  Thus,  for  higher  order  spatial  accuracy,  interpolation  stage  is  modified  without 
changing  a  Riemann  solver.  Referring  that  piecewise  constant  state  generates  first  order  spatial 
accuracy,  a  piecewise  linear  or  quadratic  distribution  is  applied  for  second  or  third  spatial 
accuracy.  From  constructing  the  adaptive  dataset,  flux  values  are  evaluated  in  the  included  cells 
in  the  adaptive  dataset.  This  research  uses  the  Finite  Volume  Method  (FVM)  to  calculate  the 
numerical  fluxes  E,  F  ;  flux  values  are  then  computed  on  all  sides  of  the  cell  with  third-order 
accuracy.  For  the  development  of  adaptive  wavelet  method  based  high  order  scheme,  we  adopt 
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TVD  and  e-MLP  scheme. 


3-4-1.  TVD  (Total  Variation  Diminishing)  limiting  process 

The  cell  interface  values  determined  by  third-order  interpolation  with  the  TVD  limiting 
function  are  then  derived  as  in  Eq.  (3-11). 


<$>l  =  <t),  +0.5max(0,min(2A®.+i,2A4>;_i,/?iA<Dj._|))  ^  (3_1  ^a) 

=  ®/+i  -0.5max(0,min(2A®  .+i,2A®.+|,^fiA0>.+|))  _  (3-1  lb) 

where  _  ®h|  ,  =  ®<+i  _  ,  a++a  =  ®i+ 2 _  ®,-+1  and  represents  a  cell 

averaged  value,  and  PL  and  PR  are 


Pl=- 


1  +  2  r1 


AO.J 

la  «  im 


3 


where  rLi  = 


Pr  ~ ' 


1  +  2  r, 


R,i+ 1 


3 


where  r} 


ra+ 1 


ao>i4  ’ 

AO  3 


(3 -12a) 
(3 -12b) 


Using  these  cell  interface  values,  spatial  discretization  is  performed  with  conventional 
schemes  such  as  AUSMPW+  or  Roe's  FDS  method. 


3-4-2.  e-MLP  (enhanced  Multi-dimensional  Limiting  Process) 

The  key  ideas  of  e-MLP  can  be  summarized  into  three  points.  First,  an  independent 
distinguishing  step,  which  is  separated  from  high  order  interpolation  and  spatial  discretization 
schemes,  is  introduced  to  a  solver.  In  this  step,  multi-dimensional  discontinuities  are  carefully 
searched  for  by  using  Gibbs  phenomena.  Subsequently,  the  computational  domain  is  divided  into 
continuous,  linear  discontinuous  and  nonlinear  discontinuous  regions  using  a  sophisticate 
method.  Second,  based  on  the  regional  information,  appropriate  limiting  criteria  are  associated 
with  each  region  in  high  order  interpolation;  in  a  continuous  region,  there  is  no  limiting  process. 
In  a  linear  discontinuous  region  such  as  a  region  of  contact  discontinuity,  the  conventional  TVD 
criterion  is  applied.  In  a  nonlinear  discontinuous  region,  e.g.,  a  region  associated  with  a  shock 
discontinuity,  MLP  is  used  for  the  removal  of  numerical  oscillation.  Third,  this  elaborate 
information  is  fed  into  a  spatial  discretization  scheme  as  well  as  a  high  order  interpolation 
scheme.  Based  on  this  information,  proper  numerical  dissipation  is  added  to  a  nonlinear 
discontinuous  region  for  robust  calculation,  as  with  the  entropy  fix  in  Roe’s  FDS. 

Through  these  ideas,  e-MLP  can  contribute  to  the  improvement  of  the  solver  in  terms  of 
accuracy,  robustness  and  efficiency.  The  brief  summary  of  3rd  order  accurate  interpolation  with 
e-MLP  is  shown  in  Eq.  (3-13-3-16). 

(a)  3rd  order  accuracy 


Pl=- 


1  +  2  r. 


'  ■>  Pr~' 


1  +  2  rB 


o  ,  —  <t> 

where  rL  ■  =  — ,+  —  r, 

+  +1 


Rj+i  <D  ** 


i+2 


i+ 1 


(3-13) 


(b)  Continuous  region:  No  limiting  function 
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=0f+O.5AAOi4> 

=  ^z+i  -  0.5^AO  ,+|  # 

(c)  Linear  discontinuous  region: 

TVD  criterion,  <fr(r)  =  max(o,min(2,2r)) 

<DZ  =  0.  +  0.5max(0,min(2A.+1,2A2  ._1,PjLAO._1)) 

l+2  1  2  1  2  5 

=  °+i  -0.5max(0,min(2A  1,2A2  !+1,PsA<D  ,))_ 

2  2  2 

(d)  Nonlinear  discontinuous  region: 

MLP  limiter,  </>(r)  =  max(0,  min(a,  ar )) 

=0),.  +0.5max(0,min(a£A  i.^AO  ,,P£A0>  ,)) 

2  2  2  ? 

=®/+i  +0.5max(0,min(«BA.+1,«BA(D.+3,PfiA<D  ,)) 

2  2  2  7 

where  aL,R  are  MLP  coefficients. 


(3 -14a) 
(3 -14b) 


(3- 15a) 
(3- 15a) 


(3- 16a) 
(3- 16b) 


After  the  interpolation  of  primitive  variables  at  the  cell  interface,  AUSMPW+  or  Roe's  FDS 
scheme  is  applied  for  the  spatial  discretization. 

3-5.  Residual  interpolation  and  time  integration 

In  cells  excluded  from  the  adaptive  dataset,  residual  interpolation  is  adopted  before  time 
integration.  As  explained  in  the  previous  section,  the  sixth-order  interpolating  polynomial  is 
needed  in  the  wavelet  decomposition  and  residual  interpolation  process  in  order  to  maintain  the 
third-order  spatial  accuracy  of  the  conventional  CFD  schemes.  Therefore,  as  in  Equation  (3-3)  in 
the  wavelet  decomposition  process,  the  sixth-order  interpolating  polynomial  is  proposed  for 
residual  interpolation,  as  shown  in  Equation  (3-17). 


~  3  25  75  75  25 

R"  =  —R ", - —R"  2  +——R’  +  — A” - - 

m  256  '*  256  ’J  128  v  128  J  256 

'J  /•)£  me-  ri  r-  /- 


— A"+4  +  —  Rn+6 

,J  256  .+4,3  256  ,+6,j 

25  n  75  n  75  „  25  n  3 

256  R,J~i  ~  256Ri  J~2  + 128  R‘J  + 128  R‘J+Z  _  256  ^  +  256  ^ 

■+U+1  =  x  (3^-4  +  3fi,V4  + 


+  3  r; 

+  77?" 


,w-21Rlj_1-21Rl1 

H-A,i  ~  27 R"-2,j  + 1747?" 


,  1  >1  A  r»  7 


(3-17) 


After  the  calculation  of  residual  values  in  the  whole  computational  domain,  time  integration  is 
performed  using  conventional  schemes. 
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Results  and  Discussion: 


1.  Compact  data  representation  using  adaptive  wavelet 

To  test  distinguish  of  region  using  wavelets,  the  test  function  is  selected  as  a  superposition  of  a 
sine  function  and  a  Gaussian  as  Eq.  (1-1). 

/ (x)  =  sin(2^x)  +  e“a(*~°-5)2  with  a  =  1 0,000  ( 1  - 1 ) 

That  is,  it  is  smooth  in  most  of  the  domain  with  a  small  interval  of  sharp  variation.  The  function 
has  the  1,025  grid  points  with  unit  interval  as  shown  in  Fig.  l-l(a).  After  5  resolution  of  the 
adaptive  wavelet  transformation,  the  remaining  number  of  grid  points  is  only  57.  Also,  as  shown 
in  Fig.  1-1  (b),  the  dataset  follows  the  local  features  of  solution.  That  is,  the  remaining  number  of 
grid  points  is  small  in  the  smooth  region  but  large  in  the  rapidly  changed  region.  This 
compressibility  of  data  is  great  advantage  of  wavelets  in  terms  of  computational  efficiency. 


(a)  Grid  points  in  an  original  dataset  (b)  Grid  points  in  an  adaptive  dataset 

Fig.  1-1.  Example  of  the  adaptive  wavelet  transformation 


Generally,  wavelet  transformation  has  procedure  of  abbreviation  of  wavelet  coefficients  that 
are  smaller  than  some  threshold  value  by  user.  Due  to  omitting  the  data,  compressed  data  and 
original  data  is  always  different,  and  this  becomes  error  of  computation.  Especially  in  CFD,  each 
numerical  scheme  has  different  order  of  accuracy  depending  on  truncation  error,  and  if  due  to  the 
compression  of  data  error  is  bigger  than  truncation  error,  it  deteriorate  the  numerical  accuracy  of 
the  conventional  schemes.  Therefore,  for  efficient  compression  of  data  with  maintaining  the 
computational  accuracy,  it  should  needed  to  research  for  thresholding  value. 

If  the  wavelet  filtering  using  the  wavelet  transformation  performed  to  CFD  data, 
computational  domain  is  divided  to  primary  and  non-primary  flow  region.  Through  the 
appropriate  flagging  grid  cells,  remained  primary  cell  is  1  and  non-primary  cell  is  0  using  tag 
after  the  wavelet  transformation.  Accordingly  original  grid  system  can  be  well  reconstructed 
using  the  only  tagged  location  as  1 .  To  reconstruction  of  grid,  simple  concept  of  one  dimensional 
grid  system  is  shown  in  Fig.  1-2. 
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Fig.  1-2  Reconstruction  of  one  dimensional  grid  system 


Data  representation  of  two  dimensional  grid  system  is  as  shown  in  Fig.  1-3.  Original  flow  data 
is  the  NACA0012  airfoil  at  that  Mach  number  is  0.8  and  angle  of  attack  is  5°.  Fig.  l-3(a)  is 
original  flow  data  and  Fig.  l-3(b)  is  represented  data  of  original  data  through  the  flagging 
wavelet  of  4th  resolution.  Fig.  l-3(c)  is  flagged  region  by  tagging  cells(l  or  0)  and  Fig.  l-3(d)  is 
error  between  original  and  represented  data.  In  this  two  dimensional  data  representation  test, 
thresholding  value  is  IE-5  and  compression  ratio  is  0.41  (compressed  data  point  /  original  data 
point).  So  if  the  threshoding  value  is  small,  error  is  small  but  compression  ratio  is  big. 


(a)  Original  pressure  distribution 


(b)  Represented  pressure  distribution 


(c)  Flagged  region  as  1  (d)  Error  of  data 

Fig.  1-3  Data  representation  using  wavelet  transformation 
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2.  Numerical  tests  of  adaptive  wavelet  method  for  CFD  solution  algorithm 
2-1.  Efficiency  of  adaptive  wavelet  in  steady  problem 

Combining  all  of  the  introduced  methods,  i.e.  the  modified  threshold  value,  the  stabilization 
technique,  residual  interpolation  and  restriction  technique,  we  checked  the  overall  enhancement 
in  computational  efficiency  of  NACA0012  flow  problems.  Here,  the  AUSMPW+  scheme  [21]  is 
used  for  spatial  discretization  and  the  LU-SGS  implicit  methods  [22]  is  for  the  time  integration. 
The  total  iteration  number,  CPU  time  gain  and  difference  in  lift  coefficient  associated  to  each 
case  are  summarized  in  table  2-1.  Figs.  2-1  show  the  pressure  contours  of  the  reference  solver 
and  adaptive  wavelet  method  with  a  2nd  level  resolution,  respectively. 


Tab 


e  2-1  Test  cases  and  results  for  the  flow  at  NACA0012  airfoil 


NACA0012 

Airfoil 

Time 

Integration 

Wavelet 

Decomposition 

Level 

Iteration 

CPU 

Time 

Time 

Ratio 

Lift 

Coefficient 

Difference 

(%) 

Subsonic 

(Ma=0.3, 

AOA=5) 

Dense  grid 

Reference 

0.5497 

LU-SGS 

(201x97) 

Reference 

3698 

213.19 

1.00 

0.5340 

2.86 

Level  1 

2046 

113.84 

1.87 

0.5331 

3.02 

Level  2 

1995 

101.77 

2.09 

0.5361 

2.47 

Level  3 

2015 

103.45 

2.06 

0.5346 

2.75 

(a)  Pressure  contours  of  the  reference  solver  (b)  Pressure  contours  of  adaptive  wavelet  method 

Fig.  2-1.  LU-SGS  time  integration  for  subsonic  flow  (Ma=0.3,  AOA=5°) 


From  table  2-1,  differences  in  lift  coefficients  have  little  relation  to  the  resolution  level  and 
change  within  a  2nd  order  of  accuracy.  Concerning  the  time  ratio,  the  adaptive  wavelet  method 
with  a  2nd  level  resolution  is  better  than  a  1st  level  while  the  difference  in  the  time  ratio  between 
the  2nd  level  and  the  3rd  level  is  slight.  There  exists  the  proper  level  according  to  the  complexity 
of  solution  and  the  number  of  grid  points  in  order  to  gain  optimum  efficiency.  If  the  number  of 
grid  points  is  not  enough  to  express  the  flow  physics  accurately,  it  is  hard  to  expect 
improvements  in  computing  efficiency  simply  by  increasing  the  resolution  level  of  an  adaptive 
wavelet  method.  In  this  case,  the  appropriate  level  is  the  second. 


19 


2-2.  Efficiency  of  adaptive  wavelet  in  unsteady  problem 

Overall  enhancement  in  computational  efficiency  of  shock-vortex  interaction  problems  [2,  23] 
is  evaluated.  Here,  the  AUSMPW+  scheme  [21]  is  used  for  spatial  discretization  and  4th  order 
Runge-Kutta  explicit  method  [16]  is  for  the  time  integration.  Fig.  2-2  shows  the  adaptive  datasets 
according  to  the  wavelet  resolution  level  at  t=15  sec.  Here,  the  dataset  follows  the  change  of  flow 
features  and  many  cells  remains  in  I (s')  near  the  vortex  and  shock  regions.  In  the  other  regions, 
the  changes  of  flow  properties  are  negligible,  and  the  remaining  cells  are  sparsely  distributed. 
The  I(s')  dataset  is  effectively  adapted  to  solution’s  features.  Also,  in  these  figures,  the 
adaptive  datasets  between  level  3  and  4  are  very  similar.  There  is  a  proper  level  of  wavelet 
resolution  according  to  the  solution  complexity.  In  this  case,  the  appropriate  level  is  4. 


Wl 

lllljp' 

iawKwJwy '  -  iuMHM 

»y* -*-g  j?  * ; 

(a)  Level  1 

(b)  Level  2 

1  :  :  ::  :  :: 

. z 

(a)  Level  3  (d)  Level  4 

Fig.  2-2  Adaptive  datasets  of  shock-vortex  interaction  problem  according  to  wavelet 

decomposition  level  at  t=15  sec 


Fig.  2-3  and  2-4  show  the  change  of  density  contours,  variation  of  the  compression  ratio,  and 
CPU  time/iteration  during  the  time  integration,  respectively.  The  CPU  time/iteration  of  the 
reference  method  was  constant  1.58  sec.  At  the  beginning  of  the  calculation,  that  of  the  new 
adaptive  wavelet  method  was  much  smaller.  However,  as  time  went  on,  the  moving  vortex 


20 


penetrated  into  and  interacted  with  the  stationary  normal  shock,  which  complicated  the  flow 
features.  Then,  the  CPU  time/iteration  of  the  new  adaptive  wavelet  method  increased  gradually 
because  more  grid  points  were  included  in  the  dataset  to  adapt  the  dataset  to  the  flow  features. 
Therefore,  the  compression  ratio  decreased  and  accompanied  the  increase  of  the  CPU 
time/iteration  of  the  adaptive  wavelet  method. 


(c)  t=11.25sec  (d)t=15sec 


Fig.  2-3  Change  of  density  contour  of  shock-vortex  interaction  problem  according  to  time 


21 


2 


- Reference  method 

—  -  Modified  w  a  ve  le  t  m  e  th  o  d 

_  1  2  3  4 

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

_kJpu^' 

1  00  200  300  400  500  600  700 

Ite  ra  tio  n 


(b)  Variation  of  CPU  time/iteration 

Fig.  2-4  Variation  of  compression  ratio  and  CPU  time/iteration  of  shock- vortex  interaction 

problem  according  to  time 


(b)  New  adaptive  wavelet  method 
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(c)  Density  contours  at  a-b  line  (y=17)  (d)  Density  contours  at  c-d  line  (y=25) 

Fig.  2-5  Density  contours  of  shock- vortex  interaction  problem  at  t=15sec;  the  wavelet 

decomposition  level  is  3. 


Overall,  the  computation  time  became  about  2.0  times  faster  at  maximum  when  the  resolution 
level  of  wavelet  was  4.  The  L2  norm  between  the  solutions  of  the  reference  solver  and  the  new 
adaptive  wavelet  method  was  2.56  x  10  7 .  The  overall  efficiency  improvement  and  the  L2  error 
according  to  wavelet  decomposition  level  are  summarized  in  Table  4.  Generally,  if  the  number  of 
cells  becomes  large,  higher  computational  efficiency  can  be  obtained  because  the  portion  of  the 
region  where  the  flow  properties  change  smoothly  increases  and  the  compression  ratio  of  the 
wavelet  method  is  enhanced.  To  confirm  this  result,  I  performed  the  same  simulations  except  for 
resizing  the  number  of  grid  points  as  241x241  and  summarize  the  results  in  Table  5.  Although 
the  order  of  L2  errors  was  very  similar  for  all  cases,  the  maximum  improvement  of 
computational  efficiency  was  reduced  from  2.0  times  to  1.74  times.  The  present  method  can 
become  more  effective  in  flow  problems  that  require  a  huge  size  of  grid  points. 

Fig.  2-5  (a)  and  (b)  show  the  density  contours  of  the  reference  and  the  new  adaptive  wavelet 
method  at  t=15sec,  respectively.  The  comparisons  of  density  contours  at  y=17(a-b  line)  and 
y=25(c-d  line)  are  presented  in  Fig.  2-5  (c)  and  (d),  respectively.  The  density  contours  are  very 
similar  to  each  other.  Also  in  Fig.  2-5  (c)  and  (d),  the  original  features  around  the  vortex  core  or 
the  shock  discontinuity  region  are  exactly  represented  by  the  adaptive  wavelet  method. 


3.  Extensive  applications  of  adaptive  wavelet  based  high  order  scheme 
3-1.  High  accurate  adaptive  wavelet  with  TVD 
3-1-1.  Shock  Sine  Wave  Interaction  Problem 

First,  the  higher  order  accurate  adaptive  wavelet  method  was  tested  with  the  one-dimensional 
shock-sine  wave  interaction  problem.  The  governing  equation  was  the  one-dimensional  Euler 
equation  and  the  initial  conditions  are  shown  in  Eq.  (3-1). 

(p,U,p)  =  (3.8571,2.6294,10.333)  where  0<x<l, 

(p,  U,  p)  =  (1  +  0.2  sin(5(x  -  5)),0,1)  where  1  <  x  <  1 0 .  (3  - 1 ) 

Here,  3rd  order  interpolation  with  TVD  limiting  function  was  utilized  for  the  evaluation  of  cell 
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interface  values.  Then,  AUSMPW+  method  was  used  for  spatial  discretization  and  3rd  order  of 
Runge-Kutta  method  was  for  time  integration,  respectively  [21,  24,  25].  And  the  number  of  grid 
points  is  set  as  mx  2Z  +1  because  of  easy  dyadic  coarsening  during  the  wavelet  transformation. 
The  value  of  8  was  set  as  l(T5 . 

In  order  to  confirm  that  the  adaptive  wavelet  method  maintains  the  numerical  accuracy  of  the 
conventional  schemes,  grid  convergence  tests  are  performed.  Fig.  3-1  shows  the  results  of  the  3rd 
order  accurate  conventional  method  and  the  adaptive  wavelet  method.  It  is  shown  in  Fig.  3  that 
the  higher  order  accurate  adaptive  wavelet  method  maintains  the  3rd  order  of  accuracy  of  a 
conventional  solver  by  setting  the  thresholding  values. 


Fig.  3-1  Results  of  grid  convergence  test. 

The  CPU  times  of  the  higher  order  adaptive  wavelet  method  are  compared  to  those  of  a 
conventional  3rd  order  accurate  CFD  solver  according  to  the  wavelet  decomposition  level  in 
Table  3-1.  In  the  case  of  higher  order  accurate  adaptive  wavelet  method  with  a  1st  level 
resolution,  the  CPU  time  rather  increases.  This  is  because  the  increase  in  CPU  time  due  to 
additional  computations  for  the  adaptive  wavelet  method  exceeds  the  decrease  in  CPU  time  from 
the  flux  evaluation  step.  On  the  other  hand,  in  the  case  of  other  resolutions,  the  decrease  in  CPU 
time  for  the  evaluation  of  fluxes  is  dominant  and  the  overall  CPU  time  decreases.  For  example, 
in  a  3rd  level  of  wavelet  decomposition  with  8001  grid  points,  the  computation  of  the  wavelet 
method  is  1.93  times  faster  than  that  of  a  conventional  solver. 


Table  3-1  Comparison  of  CPU  times  in  shock-sine  wave  interaction  problem 


T  ,  Grid 

Level  XT  u 

Number 

CPU  Time 

Original  Wavelet 

Time 

Ratio 

1001 

70.656 

73.828 

0.96 

2001 

1  4001 

139.422 

144.547 

0.97 

277.359 

286.344 

0.97 

8001 

550.453 

571.547 

0.96 

1001 

70.656 

50.531 

1.40 

7  2001 

139.422 

97.000 

1.44 

4001 

277.359 

189.594 

1.46 

8001 

550.453 

377.844 

1.46 

3  1001 

70.656 

48.547 

1.46 
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2001 

139.422 

84.891 

1.64 

4001 

277.359 

143.047 

1.94 

8001 

550.453 

285.641 

1.93 

Fig.  3-2  shows  the  solutions  of  a  conventional  solver  and  the  higher  order  adaptive  wavelet 
method  at  t=l  with  8001  grid  points.  The  original  solution  of  a  conventional  solver  and  that  of 
the  wavelet  method  are  almost  same.  In  order  to  verify  it  more  precisely,  the  density  contours  of 
a  conventional  solver  and  the  wavelet  method  with  3rd  level  resolution  are  enlarged  at  positions  A 
and  B  in  Fig.  3-3.  In  both  rapidly  changing  region  of  A  and  the  smooth  region  of  B,  the  results  of 
adaptive  wavelet  method  show  good  agreement  with  those  of  a  conventional  solver. 


Fig.  3-2  Comparison  of  density  distribution  of  the  shock-sine  wave  interaction  problem  at 

t=lsec. 


(a)  Magnified  pressure  plots  in  the  region  of  A  (b)  Magnified  pressure  plots  in  the  region  of  B 

Fig.  3-3  Detailed  density  distribution. 

3-1-2.  Shock- Vortex  Interaction  Problem 

To  assess  the  accuracy  and  efficiency  of  the  higher  order  accurate  adaptive  wavelet  method  in 
more  complicated  case,  it  is  applied  to  a  shock- vortex  interaction  problem  [23],  The  domain  is 
set  as  -20  <  x  <  5  and  -20  <  y  <  10  with  a  600x600  grid.  The  initial  velocity,  density  and  pressure 
distributions  of  a  vortex  flow  are  presented  in  Eq.  (3-2). 
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Tangential  velocity:  u0  =Mvr exp[(l  -  r2)/  2] , 

Radial  velocity: 

Vorticity:  co{r)  =  Mv(2-r2) exp[(l  -  r2 ) / 2] ,  (3-2) 

Pressure:  ^W  =  -P~Mv  exp(l-r2)f /(r_1) 
y  2 

Density:  f*r)  =]-[\-^M2v  exp(l  -  r2 . 

The  Mach  number  of  vortex  is  0.39  and  the  initial  vortex  core  is  located  at  (-5,-5).  The  normal 
shock  with  Mach  number  of  1.29  is  propagated  to  the  vortex  from  right  boundary  of  the 
computational  domain.  Fig.  3-4  shows  the  schematic  diagram  of  shock-vortex  interaction 
problem;  in  stepl,  the  shock  wave  deforms  as  the  shock- vortex  interaction  develops.  In  step2,  the 
shock  wave  deforms  more  and  more  with  the  development  of  the  interaction  and  reflected  waves 
emanate.  Finally  in  step3,  slip  lines  are  formed  and  the  lines  meet  the  reflected  wave  and  normal 
shock  at  the  triple  points.  For  the  simulation  of  this  complicated  flow  problem,  AUSMPW+ 
method  combined  with  3rd  order  interpolation  with  TVD  limiting  function  was  used  for  spatial 
discretization  and  3rd  order  of  Runge-Kutta  method  was  for  time  integration,  respectively. 


Fig.  3-4  Schematics  of  shock- vortex  interaction  problem. 
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*****  y,r.  r*>:#  \  |  * 

(c)  Level  3  (d)  Level  4 


Fig.  3-5  Adaptive  datasets  of  shock-vortex  interaction  problem  according  to  wavelet 

decomposition  level  at  t=18. 

After  the  non-dimensional  time  of  18,  the  adaptive  dataset  according  to  the  level  of  wavelet 
decomposition  are  represented  in  Fig.  3-5.  It  is  shown  in  these  figures  that  the  dataset  follows  the 
local  features  of  the  solution,  accurately;  many  cells  are  automatically  remained  near  the  vortex 
core,  normal  shock  regions,  reflected  wave  region,  and  slip  lines,  etc.  according  to  the  crucial 
flow  features.  In  the  other  smooth  regions,  the  changes  if  the  flow  properties  are  negligible  and 
the  remaining  cells  are  sparsely  distributed.  Here,  the  adaptive  datasets  between  level  3  and  4  are 
very  similar.  There  is  a  proper  level  of  wavelet  resolution  according  to  the  solution  complexity 
and  the  appropriate  level  is  3  in  this  case.  After  the  construction  of  an  adaptive  dataset,  operation 
with  high  cost  such  as  high  order  interpolation  and  flux  evaluation  are  performed  only  at  the 
remaining  cells,  which  enhance  the  computational  efficiency. 

Fig.  3-6  shows  the  density  contours  of  a  conventional  2nd  order  accurate  solver,  3rd  order 
accurate  solver  and  the  adaptive  wavelet  method  with  3  level  of  wavelet  decomposition  after  the 
non-dimensional  time  of  18.  In  this  figure,  it  is  shown  that  the  contour  of  2nd  order  accurate 
conventional  solver  is  smeared  comparing  with  that  of  3rd  order  accurate  solver.  However,  the  3rd 
order  accurate  adaptive  wavelet  method  can  maintain  the  original  features  of  a  solution;  the 
contours  of  slip  line,  reflected  wave  pattern,  etc.  are  clearly  presented  not  only  in  the  results  of 
the  3  rd  order  accurate  solver  but  also  in  the  adaptive  wavelet  method. 
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(a)  2nd  order  of  conventional  method 


(c)  3rd  order  of  adaptive  wavelet  method 
Fig.  3-6  Comparison  of  density  contours  of  shock  vortex  interaction  problem  at  t=18 


Fig.  3-7  shows  the  detailed  density  distributions  of  2nd  order  accurate  solver,  3rd  order 
accurate  solver,  the  adaptive  wavelet  method  at  y=-5  and  they  are  compared  with  the  result  of 
dense  grid  system.  In  the  result  of  2nd  order  accurate  solver,  the  vortex  and  reflected  wave  is 
dissipated  compared  with  the  3rd  order  accurate  conventional  solver  and  the  wavelet  method. 
However,  by  using  3rd  order  accurate  method,  the  density  distributions  are  clearly  represented; 
the  results  of  3rd  order  accurate  method  follows  that  of  dense  grid  system  more  accurately  though 
small  number  of  grid  points  is  used.  Therefore,  it  can  be  verified  in  this  test  that  the  adaptive 
wavelet  method  can  maintain  the  numerical  accuracy  of  the  3rd  order  accurate  conventional 
solver;  it  presents  higher  order  accurate  results  compared  with  that  of  the  2nd  order  accurate 
solver. 
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(a)  Density  distribution  at  y=-5  (b)  Detailed  density  distribution  at  A  region 


(c)  Detailed  density  distribution  at  B  region  (d)  Detailed  density  distribution  at  C  region 
Fig.  3-7  Comparison  of  density  distribution  at  y=-5. 

The  CPU  times  of  conventional  solver  and  adaptive  wavelet  method  and  L2  error  norm 
between  the  results  of  3rd  order  accurate  conventional  solver  and  adaptive  wavelet  method  are 
summarized  in  Table  3-2.  In  this  table,  the  order  of  L2  error  norm  is  very  small,  that  is,  the 
results  of  adaptive  wavelet  method  are  nearly  similar  to  that  of  the  3rd  order  accurate  solver.  On 
the  other  hand,  there  is  the  substantial  enhancement  of  the  computational  efficiency  by  using  the 
adaptive  wavelet  method;  in  case  of  the  3rd  order  accurate  solver,  the  computational  speed  is  only 
about  80%  of  the  speed  of  the  2nd  order  accurate  solver  due  to  the  operational  cost  of  high  order 
interpolation.  However,  the  computational  time  of  adaptive  wavelet  method  with  1  level  of 
wavelet  decomposition  becomes  somewhat  faster  than  the  2nd  order  accurate  solver.  The  cost  of 
high  order  interpolation  and  flux  evaluation  decreases  because  high  order  interpolation  and  flux 
evaluation  are  not  performed  at  the  excluded  cells  in  the  adaptive  dataset.  As  the  wavelet 
resolution  level  increases,  the  adaptive  wavelet  method  with  3rd  order  accuracy  becomes  more 
efficient  than  2nd  order  accurate  solver  as  well  as  3rd  order  accurate  solver;  the  computational 
time  of  adaptive  wavelet  method  is  about  1.4  times  faster  than  the  2nd  order  accurate 
conventional  solver  and  1.9  times  faster  than  the  3rd  order  accurate  solver  at  maximum  when  the 
resolution  level  of  wavelet  is  3. 

Table  3-2  Results  of  efficiency  improvements  in  the  shock- vortex  interaction  problem; 

Grid  size  is  601x601 
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Wavelet  Decomposition 

L2  Error 

CPU  Time 

Time  Ratio 

2nd  Order  Conventional 

8080 

3  rd  Order  Conventional 

10446 

0.77 

Level  1 

6.90xl0-8 

7837 

1.03 

Level  2 

8.65xl0-8 

6048 

1.34 

Level  3 

8.46xl0'8 

5603 

1.44 

Level  4 

7.42x1 0-8 

5676 

1.42 

3-1-3.  Isentropic  vortex  problem 


In  order  to  assess  the  adaptive  wavelet  analysis  in  continuous  problem,  it  was  applied  to  an 
isentropic  vortex  problem.  The  vortex  model  was  defined  as  Eq.  (3-3)  [26], 

I”  1  -  r2  " 

Tangential  velocity:  ue=Mvrex p - 


Radial  velocity:  ur  =  0 , 
Density  velocity:  p(r)  = 


(3-3) 


\  r-i 

v  2 


A 


M2exp[l-r2] 


i/(r-i) 


Pressure  distribution:  p(r)  =  pr  /  y . 


Here,  the  computational  domain  ranged  from  (-5,-5)  to  (5,5)  with  equal  grid  spacing  and  the 
vortex  core  was  located  at  (0,0).  The  vortex  Mach  number,  Mv  was  1.0  and  the  number  of  grid 
points  was  201x201 .  The  CFL  number  was  0.5  and  the  boundary  was  fixed  with  initial  values. 
For  a  simulation,  a  third-order  spatial  interpolation  with  a  MLP  limiting  function  was  utilized  for 
the  evaluation  of  cell  interface  values  and  AUSMPW+  was  for  flux  calculation.  Then  time 
integration  by  the  3rd  order  of  Runge-Kutta  method  was  applied.  After  a  non-dimensional  time  of 
50,  the  density  distributions  were  plotted. 

Through  the  thresholding,  many  cells  are  remained  near  the  vortex  core  as  shown  in  Fig.  3-8. 
Fig.  3-9  shows  the  density  distribution  of  second-order  accurate  solver,  third-order  accurate 
solver,  and  the  adaptive  wavelet  with  3  level  of  resolution  along  y=0.  The  density  distribution  of 
second  order  accurate  solver  is  damped  out  near  the  vortex  core  region.  However,  third  order 
accurate  solver  and  adaptive  wavelet  can  offer  improved  steepness  near  the  vortex  core  region. 


(a)  Level  1 
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(b)Level  2  (c)  Level  3 

Fig.  3-8  Adaptive  datasets  of  the  isentropic  vortex  flow  according  to  resolution  level. 


Fig.  3-9  Comparison  of  the  density  distributions  at  y=0. 

The  computational  time  of  second  order  solver,  third  order  solver  and  adaptive  wavelet 
method  and  the  L2  error  norm  between  the  results  of  the  third-order  solver  and  the 
multi-resolution  analysis  are  summarized  in  Table  6-3.  It  is  shown  that  the  computation  of 
multi-resolution  analysis  with  3  level  of  resolution  becomes  1.7  times  faster  than  third  order 
accurate  solver  while  the  L2  error  norms  are  nearly  small.  Hence,  multi-resolution  analysis  can 
enhances  the  accuracy  of  a  solution  compared  with  a  second  order  accurate  solver  as  well  as  it 
improves  the  computational  efficiency  of  a  third  order  accurate  solver. 

Table  3-3.  Results  of  efficiency  improvements  in  the  isentropic  vortex  problem; 


Grid  size  of201x201. 


Adaptive  wavelet 

L2  Error 

CPU  Time 

Time  Ratio 

2nd  /  wavelet(3rd) 

Second-Order  Conventional 

4081 

Third-Order  Conventional 

7054 

Level  1 

6.81  xlO"7 

5107 

1.38 

Level  2 

8.03x1  O'7 

4549 

1.55 

Level  3 

7.80x1  O'7 

4266 

1.65 
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3-2  High  accurate  adaptive  wavelet  with  e-MLP 
3-2-1  Shock- Vortex  Interaction  Problem 

First,  the  higher  order  accurate  adaptive  wavelet  method  was  tested  with  the  shock- vortex 
interaction  problem  in  Eq.  (3-2).  The  domain  is  set  as  -20<x<5  and  -20<y<10  with  a 
400  x  400  grid.  The  Mach  number  of  vortex  is  0.39  and  the  initial  vortex  core  is  located  at  (-5,-5). 
The  normal  shock  with  Mach  number  of  1.29  is  propagated  to  the  vortex  and  8  is  set  as  10"5. 
The  CPU  times  of  2nd  order  accurate  TVD  computation,  3rd  order  accurate  TVD  computation, 
3rd  order  accurate  e-MLP  computation,  and  3ld  order  accurate  e-MLP  computation  with  adaptive 
wavelet  method  are  summarized  in  Table  3-4.  In  this  table,  it  is  definitely  shown  that  3rd  order 
accurate  TVD  and  e-MLP  computations  need  more  computational  time  than  2nd  order  accurate 
TVD  computation.  Also,  e-MLP  spends  20%  of  additional  computational  cost  comparing  with 
the  result  of  3rd  order  TVD  solver  because  of  the  complexity  of  the  calculation  of  MLP 
coefficients,  aL,R  in  nonlinear  discontinuous  regions. 


Table  3-4.  Comparison  of  CPU  times  in  shock-vortex  interaction  problem. 


Level 

CPU  Time 

Time  Ratio 

2nd  order 

3rd  order 

e-MLP  (3rd) 

Wavelet 

2nd  /  wavelet(3rd) 

Level  1 

1962 

0.886 

Level  2 

1738 

2092 

2498 

1537 

1.131 

Level  3 

1478 

1.176 

However,  the  CPU  time  of  e-MLP  solver  with  1  level  of  wavelet  resolution  becomes  smaller 
than  that  of  3ld  order  TVD  solver  as  well  as  e-MLP  solver  though  it  is  still  larger  than  that  of  2nd 
order  TVD  solver.  It  is  because  the  region  which  needs  high  order  accurate  flux  evaluation 
including  e-MLP  decreases  by  the  wavelet  transformation.  The  decrease  of  CPU  time  in  flux 
evaluation  exceeds  the  increase  of  CPU  time  caused  by  wavelet  transformation  and  the 
application  of  e-MLP.  By  increasing  the  wavelet  resolution  level,  then  the  computational 
efficiency  is  substantially  enhanced  and  the  computation  becomes  faster  than  the  2nd  order 
accurate  TVD  solver;  the  CPU  time  of  wavelet  method  became  about  1.2  times  faster  than  the 
2nd  order  accurate  TVD  solver  at  maximum  when  the  resolution  level  of  wavelet  is  3. 
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■  Contnity:  5274% 

UvMir  discofitiiHity  ;  221  % 

■  Nonfciwr  discontinuity  :  253  % 

Fig.  3-10  Adaptive  dataset  by  the  wavelet  transformation  at  t=18. 

The  adaptive  dataset  by  the  wavelet  transformation  is  constructed  so  as  to  follow  these 
physical  features  of  the  solution  accurately  as  Fig.  3-10;  many  cells  are  remained  near  the  vortex 
and  shock  regions.  In  the  other  smooth  regions,  the  changes  of  the  flow  properties  are  negligible 
and  the  remaining  cells  are  sparsely  distributed.  By  the  wavelet  transformation,  the  remaining 
cells  in  the  adaptive  dataset  are  about  58%  of  the  whole  computational  domain  where  fluxes  are 
calculated.  Here,  the  remaining  regions  are  divided  into  continuous,  linear  dis-continuous  and 
nonlinear  discontinuous  regions;  53%  of  the  whole  computational  domain  is  continuous  region, 
2%  is  linear  discontinuous  region,  and  3%  is  nonlinear  discontinuous  region.  Then  in  the 
continuous  region,  there  is  no  application  of  limiting  function.  In  the  linear  discontinuous  region 
and  non-linear  discontinuous  region,  TVD  and  MLP  limiting  function  are  used  for  oscillation 
removal,  respectively.  Due  to  the  combination  of  wavelet  and  e-MLP,  the  computational 
efficiency  is  substantially  enhanced. 

Fig.  3-11  shows  the  density  contours  of  2nd  order  TVD  solver,  3rd  order  TVD  solver,  3rd  order 
e-MLP  solver  and  3rd  order  e-MLP  solver  with  the  wavelet  method  after  the  non-dimensional 
time  of  18.  In  these  figures,  it  is  shown  that  the  contour  of  2nd  order  accurate  TVD  solver  is 
smeared  comparing  with  that  of  3rd  order  accurate  solvers.  However,  the  3rd  order  e-MLP  with 
the  wavelet  method  can  maintain  the  original  features  of  a  solution;  the  contours  of  slip  lines, 
reflected  wave  pattern,  etc.  are  clearly  presented. 
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(C) 


3rd 


order  e-MLP  (d)  3  rd  order  e-MLP  +  wavelet 

Fig.  3-11  Comparison  of  the  density  distributions  at  y=0. 


Fig.  3-12  shows  the  detailed  density  distributions  of  2nd  order  TVD  solver,  3rd  order  TVD  and 
e-MLP  solver,  and  the  e-MLP  with  adaptive  wavelet  method  at  y=-5.  And  they  are  compared 
with  the  result  of  dense  grid  system  with  2nd  order  accurate  solver.  In  the  result  of  2nd  order 
TVD  computation,  the  vortex  and  reflected  wave  is  dissipated  compared  with  the  3rd  order 
accurate  solvers  and  the  wavelet  method.  However,  in  case  of  the  3rd  order  TVD,  e-MLP 
methods  and  the  wavelet  method,  the  density  distributions  follow  that  of  dense  grid  system 
accurately. 
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X 


(a)  Density  distribution  at  y=-5 


(c)  Detailed  density  distribution  at  A  region  (d)  Detailed  density  distribution  at  B  region 
Fig.  3-12  Comparison  of  density  distribution  at  y=-5. 


3-2-2.  Viscous  shock  tube  Problem 


The  second  test  is  unsteady  viscous  shock  tube  problem  [27,  28],  The  computational  domain  is 
set  to  be  0<x<l,  0<y <0.5  ,  and  the  total  number  of  grids  is  401x201.  For  boundary 
conditions,  symmetric  boundary  conditions  are  applied  at  y  =  0.5 .  Then  at  the  other  sides,  wall 
boundary  conditions  are  enforced.  The  initial  conditions  of  this  problem  are  given  as  Eq.  (3-4). 

120 

(p,u,v,p)  =  (120,0,0, - )  where  0<x<0.5, 0<y<0.5 

r 

1  2 

(p,u,v,p)  =  (1.2, 0,0, — )  where  0.5<x<l,  0<y  <0.5  (3-4) 

r 

with  y  =  1 .4,  Re  =  200,  and  p  =  constant . 


In  this  problem,  the  governing  equation  is  two-dimensional  Navier-Stokes’  equations.  Inviscid 
fluxes  are  discretized  by  AUSMPW+  &  e-MLP.  Viscous  fluxes  are  calculated  by  the  2nd  order  of 
spatial  accuracy.  For  a  time  integration,  the  4th  order  of  Runge-Kutta  method  with  CFL 
number=0.4  is  applied  and  s  is  10“5. 
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(a)  Level  1 


(b)  Level  2  (c)  Level  3 

Fig.  3-13  Compression  ratio  of  shock  tube  problem  at  each  time 

Fig.  3-13  shows  the  compression  ratio  at  each  time.  Due  to  the  moving  shock  and  interactions 
between  the  shock  and  the  boundary  layer,  the  total  compression  ratio  increase  as  time  goes 
further.  Especially,  the  compression  ratio  in  continuity  region  influences  total  compression  ratio 
most,  because  the  flow  features  become  very  complex.  However,  the  compression  ratio  in  linear 
discontinuity  region  and  non-linear  discontinuity  region  are  few  changes  of  ratio. 

Fig.  3-14  shows  the  adaptive  dataset  according  to  the  level  of  wavelet  resolution  at  t=l.  At  the 
beginning,  the  strong  shock  is  propagated  to  the  right  wall,  and  then  reflected  and  re -propagated 
again  to  the  left  accompanying  complex  shock-shear-boundary  layer  interactions.  Here,  the 
dataset  is  seen  to  follow  the  flow  feature  effectively.  Near  shock  discontinuities  and  the  vortex 
region,  there  remain  many  grid  points  to  manifest  the  crucial  features  of  solution. 
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■  Continuity  region 

■  Linear  Discontinuity  region 

■  Nonlinear  Discontinuity  region 


(a)  Level  1 


(b)  Level  2  (c)  Level  3 

Fig.  3-14  Adaptive  dataset  according  to  the  wavelet  decomposition  at  t=l. 


(c)  Adaptive  wavelet  level  3  (d)  Adaptive  wavelet  level  3 

Fig.  3-15  Density  contours  of  shock  tube  problem  at  t=l. 


In  Fig.  3-15,  the  density  contours  of  the  original  solver  and  the  adaptive  wavelet  method  are 
presented  at  t=lsec.  In  these  figures,  the  3rd  order  accurate  adaptive  wavelet  method  can  maintain 
the  original  features  of  a  solution;  the  contours  of  slip  line,  reflected  wave  pattern,  etc.  are  very 
similar  to  the  contour  of  the  3rd  order  accurate  original  solver.  Fig.  3-16  shows  the  detailed 
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density  distributions  of  3rd  order  accurate  original  solver  and  the  adaptive  wavelet  method  at 
y=0.153,  In  Fig.  3-15  (a),  the  density  distributions  of  adaptive  wavelet  method  follows  that  of 
original  solver  accurately.  But,  in  Fig.  3-15  (b),  (c),  and  (d),  the  results  of  adaptive  wavelet 
method  are  slightly  different  with  the  original  solver  because  of  thresholding  process  in  adaptive 
wavelet  algorithm. 


(a)  Density  distribution  at  y=0.077  (b)  Detailed  distribution  in  region  A 


(c)  Detailed  distribution  in  region  B  (d)  Detailed  distribution  in  region  C 

Fig.  3-16  Comparison  of  density  distribution  at  y=0.077 


The  overall  efficiency  improvement  according  to  wavelet  resolution  level  is  summarized  in 
Table  3-5.  In  this  table,  the  computational  time  with  1  level  of  wavelet  resolution  is  larger  than 
that  of  2nd  order  accurate  conventional  solver.  It  is  because  the  additional  CPU  time  of  wavelet 
transformation  is  larger  than  the  decrease  of  CPU  time  of  flux  evaluation.  By  increasing  the 
wavelet  resolution  level,  the  adaptive  wavelet  method  with  3rd  order  accuracy  becomes  more 
efficient  than  2nd  order  accurate  solver  as  well  as  3rd  order  accurate  original  solver;  the 
computational  time  of  wavelet  method  became  about  1.4  times  faster  than  the  3rd  order  of 
e-MLP  method  at  maximum  when  the  resolution  level  of  wavelet  is  3. 
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Table  3-5.  Comparison  of  CPU  times  in  shock  tube  problem;  Grid  size  (600x300) 


Total 

Comp.  Ratio 

CPU  Time  (Sec.) 

Time  Ratio 

e-MLP  (3rd) 

Adaptive  Wavelet 

Level  1 

48.21  % 

28,531 

1.181 

Level  2 

42.51  % 

33,982 

25,114 

1.353 

Level  3 

41.58% 

24,763 

1.372 

4.  Conclusions 

In  this  research,  high  order  accurate  method  for  numerical  flow  analysis  using  adaptive 
non-linear  wavelets  is  developed.  For  the  objective,  the  non-linear  wavelet  basis  is  proposed  and 
applied  to  accurate  data  representation  method.  The  CFD  solution  algorithm  using  adaptive 
wavelet  is  also  developed  and  implemented  to  high  accurate  schemes  of  TVD  and  e-MLP 
respectively.  The  proposed  adaptive  wavelet  methods  are  proved  to  be  efficient  and  preserve  the 
high  accuracy  through  the  following  numerical  exercises. 

1.  Using  the  developed  non-linear  wavelet  basis,  the  original  data  was  compressed  for  two 
applications  of  sine  function  and  a  Gaussian  and  pressure  distribution  around  the  transonic 
airfoil.  In  sine  function  and  a  Gaussian,  only  57  grid  points  can  represent  original  1025  grid 
points.  In  pressure  distribution  around  transonic  airfoil,  about  40%  of  data  can  represent 
original  data. 

2.  The  constructed  adaptive  wavelet  algorithm  is  applied  to  typical  problems  containing  the 
two-dimensional  Euler  equations.  Throughout  computations  of  a  NACA0012  airfoil  and 
shock-vortex  interaction  problems,  the  computational  efficiency  is  enhanced  as  follows.  In 
NACA0012  airfoil  problems,  by  applying  the  adaptive  wavelet  method  to  the  LU-SGS  time 
integration  method,  the  computation  is  2.1  times  faster  in  the  subsonic  case  and  1.5  times  in 
the  transonic  case  compared  with  the  results  of  the  reference  solver.  In  shock-vortex 
interaction  problem,  the  computation  becomes  2  times  faster. 

3.  The  higher  order  accurate  adaptive  wavelet  method  with  TVD  scheme  was  applied  to 
one-dimensional  shock-sine  wave  interaction,  two-dimensional  shock-vortex  interaction,  and 
the  isentropic  vortex  problems.  In  case  of  shock-sine  wave  interaction  problem,  it  was 
verified  that  the  original  numerical  accuracy  of  the  3  rd  order  conventional  solver  is  preserved 
with  enhanced  computational  efficiency.  Through  the  computation  of  shock-vortex 
interaction  problem,  it  was  shown  that  adaptive  datasets  are  constructed  so  as  to  capture  local 
features  of  solution  automatically.  Also,  the  higher  order  accurate  adaptive  wavelet  method 
speeds  up  the  simulation  1 .4  times  faster  than  that  of  the  2nd  order  accurate  solver  while  the 
3rd  order  numerical  accuracy  of  a  solution  is  maintained.  Also,  the  proposed  3rd  order 
method  was  about  1.65  times  faster  than  2nd  order  conventional  solver  while  preserving  3rd 
order  of  spatial  accuracy  in  isentropic  vortex  problem. 

4.  The  higher  order  accurate  adaptive  wavelet  method  with  e-MLP  scheme  was  applied  to 
two-dimensional  shock-vortex  interaction  and  the  viscous  shock  tube  problems.  In  the 
shock-vortex  interaction  problem,  the  proposed  3rd  order  adaptive  wavelet  method  can 
present  much  better  efficiency  (1.2  times  faster)  and  accuracy  compared  with  2nd  order 
accurate  computation  by  decreasing  the  computational  time  of  3rd  order  accurate  e-MLP 
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method.  In  the  viscous  shock  tube  problem,  the  adaptive  wavelet  method  with  3rd  order 
accuracy  becomes  more  efficient  than  2nd  order  accurate  solver  as  well  as  3rd  order  accurate 
original  solver;  the  computational  time  of  wavelet  method  became  about  1 .4  times  faster  than 
the  3  rd  order  of  e-MLP  method. 

The  developed  high  accurate  adaptive  wavelet  method  was  confirmed  to  improve  the 
computational  efficiency  of  the  reference  solver  and  to  improve  the  numerical  accuracies 
compared  to  the  similar  computation  workload  problem.  These  capabilities  of  the  proposed 
algorithm  were  successfully  enhanced  by  adopting  three  techniques  for  the  adaptive  wavelet 
method. 

1 .  The  threshold  value  is  modified  to  consider  the  grid  spacing  and  the  size  of  temporal  step.  As 
a  result,  it  is  realized  that  the  adaptive  wavelet  method  preserves  the  spatial  and  the  temporal 
accuracies  of  the  reference  solver. 

2.  The  stabilization  technique  is  applied  to  the  construction  process  of  an  adaptive  dataset  in 
order  to  improve  the  compression  ratio.  Also,  residual  interpolation  is  performed  at  the  11 
time  step,  not  at  the  n  + 1  time  step.  It  enables  the  process  of  arbitrary  addition  of  adjacent 
cells  to  a  I(s')n  dataset  to  be  excluded  from  the  wavelet  process.  Therefore,  the  I  (s')" 
dataset  can  be  automatically  constructed  without  the  user’s  arbitrary  adjustment  of  the 
threshold  value. 

3.  In  the  time  integration,  if  the  variations  of  the  flow  variables  are  smaller  than  the  threshold 
value  s' ,  they  are  controlled  by  adopting  a  weighting  factor.  This  restriction  technique 
enhances  the  convergence  rate  of  steady  state  calculations. 
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Efficiency  Enhancement  in  High  Order 
Accurate  Euler  Computation  Via  AWM 


Hyungmin  Kang.  Dongho  Lee,  Dohyuog  Lee.  Dochan  Kwak  and  John  Seo 


I  Introduction 

Recently,  in  various  aerospace  transportation  and  exploration  system,  l here  have 
been  strong  demands  for  accurate  compulations  of  flow  features  including  super¬ 
sonic  or  hypersonic  flow  phenomena.  These  demands  for  accurate  computations 
oblige  the  developments  in  Computational  Fluid  Dynamics  1CFD)  modeling  and 
simulation  technologies.  As  a  result.  CFD  tools  have  been  widely  used  in  the 
aerospace  engineering  and  have  become  indispensable  in  the  part  of  design  stage 
of  aerospace  vehicles. 

However,  the  current  level  of  technology  does  not  meet  the  requirement  of  the 
CFD  field  yet.  For  the  enhancement  of  the  prediction  capability  of  CFD  tools,  it  is 
essential  to  perform  high  fidelity  simulations  with  many  grid  points  and  high  order 
accurate  algorithms.  Then,  the  increment  of  computational  cost  is  inevitable,  espe¬ 
cially  in  unsteady  calculations.  In  order  to  extend  the  availability  of  high  fidelity 
simulations,  it  is  necessary  to  alleviate  the  computational  cost  through  the  advance¬ 
ment  of  CFD  algorithms. 

From  the  viewpoint  of  computational  efficiency,  the  adaptive  wavelet  method  can 
be  a  good  remedy.  The  crucial  idea  of  the  adaptive  wavelet  method  is  as  follows:  in  a 
CFD  dataset,  the  majority  of  region  is  smooth  except  some  rapidly  changing  region 
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such  as  shock,  boundary  layer,  vortices,  etc.  Then  the  wavelet  coefficients  computed 
from  decomposition  process  are  compared  with  the  threshold  value,  £*.  If  these  val¬ 
ues  are  larger  than  c',  points  become  remained  in  an  adaptive  dataset.  Through  these 
processes,  tlie  adaptive  wavelet  method  can  present  the  automatic  adaptation  of  a 
dataset  based  on  the  local  feature  of  a  CFD  solution  [1-3].  Then,  flux  evaluation  is 
performed  only  at  the  sparsely  remaining  points  in  an  adaptive  dataset,  which  sub¬ 
stantially  reduces  the  computational  time.  However,  by  the  thresholding  process, 
additional  errors  are  added  to  a  solution  and  may  cause  the  deterioration  of  the  nu¬ 
merical  accuracy  of  conventional  CFD  schemes  if  the  additional  errors  are  larger 
than  the  truncation  errors  of  the  spatial  discretization  and  time  integration  methods. 
To  overcome  this  defect,  Kang  ei  a 1.  present  die  modified  threshold  value  which 
considers  the  temporal  accuracy  as  well  as  die  spatial  accuracy  of  the  conventional 
CFD  schemes.  And  they  applied  it  to  the  2nd  order  accurate  unsteady  flow  problems 
[2,3]. 

In  diis  paper,  the  objective  is  the  high  order  accurate  Euler  computation  with 
enhancing  the  computational  efficiency.  For  this  purpose,  3rd  order  accurate  Euler 
solver  is  constructed  with  AUSMPW+  spatial  discretization  method  and  3rd  order 
Runge-Kutta  lime  integration  scheme  [4,5],  And  Mu  Hi -dimensional  Limiting  Pro¬ 
cess  (MLP)  is  used  in  order  to  remove  numerical  oscillation  [6].  This  3rd  order  ac¬ 
curate  Euler  solver  is  combined  with  the  adaptive  wavelet  method  as  follows:  first, 
die  threshold  value  is  modified  to  be  adequate  for  maintaining  3rd  order  spatial  and 
temporal  accuracy  of  a  current  CFD  solver.  Second,  die  general  adaptive  wavelet 
transformation  procedure  is  changed  by  adopting  residual  interpolation.  Through¬ 
out  these  processes,  the  higher  order  accuracy  of  a  conventional  solver  is  conserved 
and  tlie  computational  cost  is  substantially  reduced.  In  order  to  demonstrate  the 
efficiency  and  tlie  accuracy  of  the  developed  method,  tlie  method  is  applied  to  a 
complicated  shock-vortex  interaction  problem  [7]. 


2  Numerical  Algorithm  of  the  Adaptive  Wavelet  Method 


In  this  research,  the  twro-dimensional  Euler  equations  are  used  as  the  governing 
equations  of  unsteady  flowr  problems.  The  generalized  coordinate  transformed  two- 
dimensional  Euler  equations  are  w  ritten  as  Eq.  { 1 ). 


dQ 

dt 

with  Q=^.E=]j[$,Q 


drj1-  R‘-j’ 

+  §yF],  F  =  j[q,Q  +  1]XE  +  qyF\. 


(I) 


For  constructing  3rd  order  accurate  conventional  solver,  AUSMPW+  with  MLP  lim¬ 
iting  function  is  used  for  spatial  discretization  and  3rd  order  Runge-Kutta  method  is 
for  time  integration  [4-6 J. 
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Then,  by  the  decomposition  process,  the  estimation  of  flow  variables  is  per¬ 
formed.  If  we  assume  that  (/,/)  cell  is  even  numbered  cell  and  (1+  1  J)7  (ij  + 1} 
and  (/  +  1 ,  j  +  1 )  are  odd  numbered  cells,  the  values  at  even  numbered  cell  is  saved 
to  a  point  in  the  coarser  level  grid.  At  odd  numbered  cells,  we  approximate  the  orig¬ 
inal  values  by  interpolating  polynomial.  In  this  research,  6th  order  of  interpolating 
polynomial  is  used  for  maintaining  the  3rd  order  accuracy  of  a  conventional  solver. 
The  equations  with  6th  order  of  accuracy  for  the  approximation  procedure  of  two 
dimensional  flow  problems  are  presented  in  Eq.  (2). 


+3ST-V+2  -  VQU.H2  + 1 74 e?J+2  +  -  71QU.h2  +  3 £ 

+20-2,j+4-27C7J+4-27f2V2  J+4  +  2g?+4J+4  +  3Q?j+6+3G?+2 

1/+6 


Then,  difference  values  between  original  values  and  approximation  values  are 
calculated  and  compared  with  threshold  value.  In  this  research,  in  order  to  maintain 
the  accuracy  of  a  conventional  solver,  the  modified  threshold  value,  ey  is  defined  as 


Eq.  (3). 


ef  =  min\e ,  max{  (Ax)3 , max(  \  U j  f  | V\)  ■  CFL '?  -  (Ax)3)] . 


(3) 


If  difference  value  is  larger  than  the  point  is  included  in  the  adaptive  dataset; 
if  not,  the  point  is  excluded  from  the  dataset.  Throughout  this  process,  the  dataset 
is  adapted  to  the  flow  features  while  maintaining  the  numerical  accuracy  of  con¬ 
ventional  schemes.  After  thresholding  process,  flux  values  are  only  calculated  at  the 
included  cells  in  the  adaptive  dataset  by  AUSMPW+  with  MLP  limiting  function 
with  3rd  order  accuracy.  At  excluded  cells  in  the  dataset,  residual  values  are  inter¬ 
polated  by  using  the  same  interpolating  polynomial  as  Eq.  (2).  And  3rd  order  of 
Runge-Kutta  time  integration  is  performed  on  the  whole  computational  domain. 


3  Numerical  Test  and  Discussion 

In  order  to  assess  the  numerical  accuracy  and  efficiency  of  the  higher  order  accurate 
adaptive  wavelet  method,  it  is  applied  to  a  shock-vortex  interaction  problem  [7].  The 
computational  domain  is  set  as  —20  <  x  <  5  and  —20  <y  <  10  with  a  600  x  600 
grid.  The  Mach  number  of  vortex  is  0.39  and  the  initial  vortex  core  is  located  at  (- 
5,-5).  This  vortex  is  propagated  to  a  stationary  normal  shock  with  the  Mach  number 
of  1.29.  Here,  e  is  set  as  1G-5. 
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(a)  Conventional  method 


(b)  Modified  wavelet  method 


Fig.  1  Bensily  contours  of  shock-vortex  interaction  problem  at  t=ISsec:  the  wavelet  decomposi¬ 
tion  level  Is  3. 


Figure  1  shows  the  density  contours  of  a  conventional  3rd  order  accurate  solver 
and  the  adaptive  wavelet  scheme  after  the  non-dimensional  time  of  18.  It  is  shown 
that  they  are  almost  same.  The  adaptive  wavelet  method  can  maintain  the  original 
features  of  a  solution  such  as  the  vortex  pattern,  the  position  and  strength  of  shock, 
etc.  Figure  2  show  s  the  detailed  density  distributions  of  2nd  order  accurate  solver. 
3rd  order  accurate  solver,  the  adaptive  wavelet  method  at  y=-5  and  they  are  com¬ 
pared  with  the  result  of  dense  grid  system.  In  case  of  2nd  order  accurate  case,  the 
distribution  is  smeared  compared  with  3rd  order  accurate  cases.  However,  the  re¬ 
sults  of  3rd  accurate  solver  and  3rd  order  accurate  adaptive  wavelet  method  present 
more  accurate  distributions  than  2nd  order  accurate  solver. 

Figure  3  represents  the  adaptive  dataset  w  ith  w  avelet  decomposition  level=3.  The 
dataset  follow  s  the  flowr  features  accurately  and  many  cells  are  remained  near  the 
vortex  and  shock  regions.  In  tlie  other  smooth  regions,  the  changes  of  the  flowr  prop¬ 
erties  are  negligible  and  the  remaining  cells  are  sparsely  distributed.  Fluxes  are  cal¬ 
culated  only  at  the  remaining  points,  which  enhances  the  computational  efficiency 
as  shown  in  Table  I .  In  this  table,  the  3rd  order  accurate  adaptive  wavelet  method 
presents  better  computational  efficiency  than  the  2nd  order  accurate  solver.  There¬ 
fore,  the  present  adaptive  wravelet  method  can  present  high  order  accurate  results 
w  ith  better  efficiency  than  conventional  CFD  solvers. 


4  Conclusion 

Throughout  this  research,  the  higher  order  accurate  adaptive  wavelet  method  is  pro¬ 
posed  for  the  flow  compulations  which  need  higher  order  accuracy  and  efficiency. 
For  this  purpose,  6th  order  accurate  interpolating  polynomial  is  implemented  to  the 
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(a)  Comparison  of  density  plot  at  y=-5 


Fig.  2  Comparison  of  density  plots  of  shock-vortex  interaction  problem  at  t=18sec;  the  wavelet 
decomposition  level  is  3. 


wavelet  decomposition  process.  Also,  die  threshold  value  is  modified  in  order  to 


Table  t  Results  of  efficiency  improvements  and  L2error  for  the  shock- vortex  interaction  problem 


Wavelet  Level 

CPU  Time 

Time  Ratio 

Third  Order 

1241  L3 

Level  1 

9086.0 

1.37 

Level  2 

6798.5 

1.83 

Level  3 

6020.1 

2.06 

Level  4 

6053.6 

2.05 

Second  Order 

8079.2 

1.54 
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Fig.  3  Adaptive  datasets 
of  shock -vortex  interaction 
problem  with  wavelet  decom¬ 
position  Ieve1=3  alt=t8 


maintain  the  3rd  order  accuracy  of  conventional  CFD  schemes.  This  higher  order  ac¬ 
curate  adaptive  wavelet  method  is  applied  to  the  shock -vortex  interaction  problem. 
And  it  can  present  much  better  computational  efficiency  than  a  2nd  order  accurate 
scheme  while  maintaining  the  higher  order  numerical  accuracy. 
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Abstract 

An  adaptive  wavelet  transformation  method  with  high  order 
accuracy  is  proposed  to  enable  efficient  and  accurate  flow 
computations.  The  mam  algorithm  is  as  follows;  3rd  order  of 
wavelet  decomposition  and  thresholding  are  applied  to  a 
CFD  dataset  and  the  positions  of  crucial  features  in  the 
dataset  are  searched  with  maintaining  the  original  numerical 
accuracy  of  a  conventional  solver  and  the  dataset  is  auto¬ 
matically  adapted  to  local  features  of  a  solution.  After  the 
wavelet  transformation.  3rd  order  spatial  and  temporal 
accurate  high  order  interpolation  scheme  with  Multi¬ 
dimensional  Limiting  Process  (MLP)  are  performed  only  at 
the  included  points  m  the  adapted  dataset,  hi  the  other 
points,  high  order  of  mterpolation  method  is  utilized  m 
order  to  construct  residual  values.  Hus  high  order  mter¬ 
polation  scheme  with  high  order  adaptive  wavelet  trans¬ 
formation  was  applied  to  unsteady  Euler  flowr  computations. 
Through  these  processes,  both  computational  efficiency  and 
numerical  accuracy  are  guaranteed  m  case  of  high  order 
accurate  unsteady  flow'  computations. 

1  Introduction 

Recently,  there  have  been  strong  demands  for  accurate  flow' 
computations  including  supersonic-hypersonic  phenomena 
based  on  the  enhancement  of  analysis  systems.  Especially 
m  various  aerospace  transportation  and  exploration  systems, 
the  prediction  capability  of  flow'  phenomena  becomes  im¬ 
portant  for  cost  reduction  during  die  part  of  design  stage  of 
aerospace  vehicles.  For  the  satisfaction  of  these  demands, 
accurate  modeling  algorithms  and  simulation  technologies 
with  vast  size  of  gnds  have  been  forced  in  Computational 
Fluid  Dynamics  (CFD).  As  a  result.  CFD  algorithms  have 


been  widely  used  m  various  filed  including  die  aerospace 
engineering. 

However,  the  current  level  of  technology'  is  often  in¬ 
sufficient  for  the  satisfaction  of  the  requirement  of  the  CFD 
field.  Still,  complexity'  of  the  flow'  phenomena  and  shape  of 
die  models  becomes  the  obstruction  m  terms  of  accuracy 
and  efficiency.  To  enhance  the  prediction  capability  of  CFD 
algorithms,  high  fidelity  algorithms  with  many  grid  points 
are  essential;  the  increment  of  computing  cost  is  inevitable, 
especially  m  unsteady  calculations.  Hence,  it  is  necessary’  to 
accurately  calculate  die  flowr  phenomena  and  efficiently 
alleviate  die  computational  cost  through  die  advancement  of 
the  CFD  algorithms. 

From  the  \iewpomt  of  computational  efficiency,  the  adap¬ 
tive  wavelet  method  can  be  a  good  remedy’  [1-3]. The  crucial 
idea  of  the  adaptive  wavelet  method  is  as  follow's;  generally, 
the  usage  of  dense  grid  system  m  the  whole  computational 
domain  becomes  a  waste  of  resources  because  the  majority 
of  the  domain  consists  of  smooth  regions.  A  dense  grid 
system  is  only  needed  in  the  rapidly  changing  regions  such 
as  shock  wraves,  boundary  layers,  etc.  Through  the  wravelet 
transformation,  wavelet  coefficients  become  larger  than 
threshold  value  s  m  these  crucial  regions.  Then,  the  dataset 
is  automatically  adapted  to  the  local  features  of  a  solution 
[1-3].  Flux  values  are  calculated  only’  at  the  remaining 
pomts  m  an  adaptive  dataset  wluch  reduces  the 
computational  tmie. 

However,  the  accuracy  of  the  adaptive  wavelet  method  is 
sensitive  to  s  value.  By  the  wavelet  transformation,  addi¬ 
tional  numerical  errors  might  be  added  to  the  solution  and 
the  numerical  accuracy  of  a  conventional  CFD  solver 
becomes  deteriorated  if  the  order  of  additional  error  is 
larger  than  that  of  truncation  errors  of  the  spatial  discrete- 
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zation  and  time  integration  method.  In  order  to  resolve  this 
defect,  modified  threshold  value  is  presented  by  considering 
the  spatial  and  temp  oral  accuracies  of  numerical  schemes  [2, 

3]- 

In  this  paper,  our  main  objective  is  to  efficiently  compute 
the  Euler  flows  with  higher  order  of  accuracy.  For  this 
purpose.,  higher  order  of  wavelet  trail  sfonnation  process 
applicable  to  the  3ld  order  accurate  solver  is  constructed  as 
follows;  first  the  6th  order  of  interpolating  polynomial  is 
introduced  to  wavelet  decomposition  process.  Second,  the 
threshold  value  is  modified  in  order  to  maintain  the  3rd 
ci  der  of  spatial  and  temporal  accuracies  of  a  solver.  This 
adaptive  wavelet  method  is  fed  into  the  3  rd  order  accurate 
Euler  solver  and  applied  to  several  numerical  tests.  Through 
the  application  of  the  higher  order  accurate  adaptive 
wavelet  method,  the  accuracy  of  the  conventional  solver  is 
conserved  and  die  cost  is  substantially  reduced. 

This  paper  is  organized  as  follows;  after  the  introduction, 
the  overall  adaptive  wavelet  method  including  the  higher 
order  of  decomposition  and  thresholding  is  described.  Then, 
to  assess  the  developed  method,  several  numerical  ex- 
penments  are  performed  such  as  one-dunensional  shock- 
sine  interaction  pr  oblem  and  two-dimensional  shock-vortex 
inter  action  problem.  Finally,  a  conclusion  is  drawn. 


2  The  overall  implementation  of  the  adaptive 
wavelet  method 


Hie  two  dimensional  Euler  equation  is  written  as  follows: 

=o,  CD 


dO  dE  dF 
—  +  —  +  — - 


dt  dy 


with  Q  - 
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where  all  properties  and  governing  equations  are  non- 
dmiensionahzed.  By  the  generalized  coordinate  trans¬ 
formation,  Eq.  (2)  is  derived. 


dO  rdE  dF 

— ~  —  — [ - 1 - ] 

dr  dij 


-Ks* 


(2) 


with  Q  =  %  ,  E  =  +  &E  +  £0]  311(1 

F =^j\JhQ+TlxE+lyF\ 

Hie  fiowr  chart  m  figure  1  shows  the  overall  processes  of 
the  implementation  of  the  adaptive  wavelet  method  to  a 
two-dimensional  Euler  solver. 


Figure  1 :  Overall  procedure  of  the  adaptive  wavelet  method. 
2.1  Wavelet  transformation 

At  first,  assume  a  two-dimensional  dyadic  gridset  as  figure 
2.  Here,  positions  of  the  symbol  Q  are  even  numbered 
gilds  and  the  other  positions  of  symbols  A ,  and  X  are 
odd  numbered  gnd  pomts.  Hien.  values  at  even  numbered 
cell  are  saved  to  a  pomt  in  the  coarser  level  grid.  At  odd 
numbered  cells,  6th  order  of  mteipolatmg  polynomial  is 
used  for  maintaining  the  3rd  order  accuracy  of  a  con¬ 
ventional  solver  as  shown  in  Eq.  (3). 


□  : 


0+i  j  = 


^(3^-250^  +  1500", 
+  1500"  2j—  25gfH  >+  3£SUj) 


Q '  U  =  —  (30% 4  —  250 j-2  +1500, 
+ 1 5O0^j+2  —  250"j-+4  +  30y+6> 


(3) 


Qhj  = 

-  270T+3J-2  +  20%, -2  +  30"-4  J -270%, 

+ 1 740,  + 1 740%, -  270+4j  +  30%, 

+  30%,, +2  -  270%, +2  + 1740", -+2  + 1740"  2  J+2 

-  270-4  ,+j  +  30+t2,-2  +  20-2.J+4  -  270, +4 

-  270-2 J+4  +  20-4, +4  +  30, +6  +  30+2j+% 
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Figure  2:  Example  of  two-dimensional  dyadic  grid. 

The  difference  values  are  calculated  as  Eq.  (4). 

□:<w=sr+u-a"u> 

A:  d"j+i  =  Qllj+i  - Qjj+i ,  (4) 

X  -  ^,-lJ+l  =  2"+L/+l  ~  0N-1J+1 


lmiitmg  function  may  damage  the  solution  s  accuracy  m  the 
contmuous  region  though  it  is  mdispensible  for  the  os¬ 
cillation  removal  due  to  higher  order  accurate  scheme  m  the 
discontinuous  region.  Also,  it  is  a  waste  of  the  computation¬ 
al  resources  to  use  limiting  function  m  contmuous  region. 

In  order  to  resolve  this  problem,  the  adaptive  dataset  is 
divided  into  two  regions:  contmuous  and  discontinuous 
regions  [7].  In  the  contmuous  region,  there  is  no  application 
of  limiting  process  and  only  higher  order  of  interpolation  is 
applied.  On  the  other  hand,  m  discontinuous  region,  MLP  is 
utilized  for  the  removal  of  numerical  oscillation.  Through 
this  switching  of  limiting  function,  accurate  and  efficient 
calculation  of  flux  values  is  possible. 

2.3  Residual  interpolation 

At  excluded  cells  m  the  dataset,  residual  values  are 
interpolated  by  using  the  6*  order  of  interpolating 
polynomial  as  Eq.  (6). 

□  &1J  =  ^(3^,;  -25C2j.  +  150^; 

+  150*"+2J  -  25R”+ij  +  3  *,”+6j) 


This  wavelet  decomposition  process  is  performed  by  multi- 
resolution  and  until  the  coarser  level  of  gridset  is  obtamed. 

After  wavelet  decomposition,  the  wavelet  coefficients  are 
compared  with  the  threshold  value;  if  the  values  are  larger 
than  the  threshold  value,  the  pomts  become  remained  m  an 
adaptive  dataset.  Through  these  processes,  the  dataset 
follow’s  the  local  feature  of  a  CFD  solution.  However,  the 
thresholding  process  inevitably  results  in  the  additional 
0(s)  numerical  errors  because  of  the  data  loss  below’  O(s) . 
If  O(s)  error  is  larger  than  the  truncation  errors  of  the 
spatial  discretization  and  time  integration  method,  the 
accuracy  of  a  conventional  CFD  scheme  should  be  de¬ 
teriorated.  In  order  to  overcome  this  defect,  the  thresholding 
method  is  modified  by  considering  the  order  of  truncation 
errors  as  Eq.  (5). 

s  -  min[£\  max( Ax3 ,  CFZ?  Ax3 )] .  (5) 

If  a  wavelet  coefficient  is  smaller  than  s'  at  some  pomt, 
the  grid  pomt  doesn’t  belong  to  the  dataset.  In  the  other 
case,  the  pomt  is  included  m  the  dataset;  an  adaptive 
dataset  is  constructed  according  to  the  local  features  of  the 
dataset  with  3rd  order  of  accuracy.  Through  these  higher 
order  of  wavelet  decomposition  and  modified  threshold 
value,  3ri  order  of  spatial  and  temporal  accuracies  of  a 
conventional  CFD  solver  are  guaranteed. 

2.2  Flux  evaluation 

After  the  wavelet  transformation,  flux  values  are  only 
calculated  at  the  mcluded  cells  m  the  adaptive  dataset  with 
3rd  order  of  accuracy.  In  order  to  construct  a  3rd  order 
accurate  solver,  AUSMPW+  scheme  is  utilized  [4,  5].  Also, 
for  oscillation  removal  due  to  numencal  discontinuities, 
Multi-dmiensional  Limiting  Process  (MLP)  is  applied  [6]. 
However,  there  is  no  need  of  the  application  of  limiting 
flmction  to  the  contmuous  regions;  dissipation  due  to 


+150^2-25^+3^) 

&l,j  =^(3^",-4  +  3*,>2j-4  +  Uf-Xj-l-nKj-l 

~  ^Kij-2  +  23&4J-2  +  3*,"-4.;  -  27*".2.> 
x .  +17<;  +174*,",,.;  -  27*" 4j  +  3  Jfoj 

+  3/£4.;+2  -  27*?_2,,+2  + 174*^  + 174*?+2j+2 

-  27/£  4.j+2  +  iKij-2  +  2*,"-2j+4  -  277^,4 

-  21 2  ;+4  +  2J2&4J+4  +  3  R? J+6  +  37^+2  j+6) 

By  using  this  6th  order  of  interpolating  polynomial,  the  error 
caused  by  interpolation  becomes  smaller  than  the  truncation 
errors  of  the  numencal  schemes.  After  constructing  residual 
distributions  m  die  whole  computational  domain.  3rd  older 
of  Rimge-Kutta  tmie  integration  is  performed  [4,  5]. 

3  Numerical  results  and  discussion 

In  order  to  assess  die  accuracy  of  the  higher  order  accurate 
wavelet  transfonnation,  it  is  applied  to  shock-sine  wave 
interaction  problem.  The  initial  conditions  for  the  shock- 
sine  wave  interaction  problem  are  given  m  the  following  Eq. 
(7).  Here,  s  is  set  as  10*5. 

(p,U,p)  =  (3.8571,2.6294,10.333)  where  0<x<l, 

(/?, U,p)  =  (1  +  0.2sin(5(x - 5)), 0.1)  where  1  <  x  <  10  .(7) 

The  CPU  times  of  the  higher  order  adaptive  wavelet  method 
are  compared  to  tiiose  of  a  conventional  3rd  order  accurate 
CFD  solver  according  to  the  resolution  level  m  Table  1. 
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Table.  1  Comparison  of  CPU  tunes  in  skock-sine  wave 
interaction  problem. 


Ltvd 

Grid 

CPU  Time 

Compression 

Number 

Original 

Wavelet 

Ratio 

1001 

70.656 

73.S2B 

1.672 

2001 

130,422 

144.547 

1.756 

1 

4001 

277,350 

286,344 

1,775 

8001 

550.453 

571.547 

i,7B9 

1001 

70.856 

50,531 

2,304 

2001 

139.422 

97.000 

2.54B 

2 

4001 

277,350 

189,504 

2,649 

B001 

550.453 

377,844 

2,695 

1001 

70,656 

48.547 

2.3B1 

2001 

130.422 

84,001 

2.829 

4001 

277.359 

143.047 

3,396 

0001 

550,453 

285.641 

3,507 

1001 

2001 

70.656 

139.422 

83,984 

2,857 

4001 

277.350 

140,504 

3,413 

8001 

550.453 

272.B13 

3,675 

In  the  case  of  the  adaptive  wavelet  method  with  a  1st  level 
resolution,  the  CPU  tune  rather  increases.  This  is  because 
the  increase  m  CPU  time  due  to  additional  computations  for 
the  adaptive  wavelet  method  exceeds  the  decrease  in  CPU 
time  from  the  flux  evaluation  step.  On  the  other  hand,  in  the 
case  of  the  adaptive  wavelet  method  with  a  2  nd  level 
resolution,  the  decrease  in  CPU  tune  for  the  evaluation  of 
fluxes  is  dominant  and  the  overall  CPU  tmie  decreases.  For 
example,  in  a  2nd  level  resolution  with  8001  grid  points,  the 
computation  of  the  wavelet  method  is  1 .46  times  faster  than 
that  of  a  conventional  solver. 

Figure  3  shows  the  solutions  of  a  conventional  solver  and 
the  higher  order  adaptive  wavelet  method  at  t=l.  The 
original  solution  of  a  conventional  solver  and  that  of  the 
wavelet  method  are  almost  same.  In  order  to  verify  it  more 
precisely,  the  density  contours  of  a  conventional  solver  and 
the  wavelet  method  are  enlarged  at  positions  A  and  B  and 
compared  in  figure  4.  At  A  and  B  regions,  the  results  of 
adaptive  wavelet  method  show  good  agreement  with  those 
of  a  conventional  solver. 


Figure  3:  Density  distribution  of  the  shock-sine  wave 
interaction  problem  at  t=l . 


(b)  B  region. 

Figure  4:  Detailed  density  comparison. 

Finally,  in  older  to  confirm  the  spatial  accuracy  of  this 
method,  grid  convergence  tests  are  performed  and  the 
results  are  shown  in  figure  5.  From  figure  5,  the  higher 
order  accurate  adaptive  wavelet  method  maintains  the  3]d 
order  of  accuracy  of  a  conventional  solver  by  setting  the 
thresholding  values  as  Eq.  (5). 


Figure  5:  The  result  of  grid  convergence  test. 

To  assess  the  accuracy  and  efficiency  of  the  higher  order 
accurate  adaptive  wavelet  method  in  more  complicated  case, 
it  is  applied  to  a  shock-vortex  interaction  problem  [8]  The 
domain  is  set  as  —20  <  x  <  5  and  -20  <  y  <  10  with  a 
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400X400  grid  The  initial  velocity,  density  and  pressure 
distributions  of  a  vortex  flow  are  presented  in  Eq.  (8)  [8]. 

T angenf i al  velocity:  ue  =  Mj r  ©cp  [( 1  -  r 2)  /  2] , 

Radial  velocity:  ur  =  0  , 

Voriidty:  tfXV)  =  My  (2-r1)  exp[(l  —  r 3  )  /  2] ,  (8) 

Pressure:  p(r)  = “ [1  -  — — -Mv3  exp(l  -  r2)]r/tr_1) , 

7  2 

Density:  p(r)  =  —[1  -  — — -  My  exp(l  -  r 2  )]!  ^  _1? . 

7  2 

Hie  Mach  number  of  vortex  is  0.39  and  the  initial  vortex 
core  is  located  at  (-5,-5).  The  normal  shock  with  Mach 
number  of  1.29  is  propagated  to  the  vortex  and  £  is  set  as 
10'5. 

The  overall  efficiency  improvement  according  to  wavelet 
resolution  level  is  summarized  in  Table  2.  In  this  table,  the 
computational  time  with  1  level  of  resolution  is  larger  than 
that  of  2Ed  order  accurate  conventional  solver.  It  is  because 
the  additional  CPU  time  of  wavelet  transformation  is  larger 
than  the  decrease  of  CPU  time  of  flux  evaluation,  as 
mentioned  in  Hie  previous  example.  By  increasing  the 
wavelet  resolution  level  the  adaptive  wavelet  method  with 
3rd  order  accuracy  becomes  more  efficient  than  2nd  order 
accurate  solver  as  well  as  3rd  order  accurate  solver;  the 
computational  time  became  about  2.0  tunes  faster  at 
maximum  when  the  resolution  level  of  wavelet  is  3. 

Table.  2  Comparison  of  CPU  tunes  m  shock-vortex 
interaction  problem. 


level 

CPU  Time 

Computational 

Efficiency 

[OTtm jGm'  HU  hMTH.Wai  l 

Conventional 
(Z"1  Qider) 

Conventional 
(3ri  order} 

Wavelet 
(3nJ  order) 

1 

3079,2 

12411.3 

90S6.D 

1.37 

2 

6798.5 

1,63 

3 

6020.1 

2,06 

4 

6053.6 

2.05 

Figure  6  represents  the  adaptive  dataset  with  wavelet 
decomposition  level =3 .  The  dataset  follows  the  flow 
features  accurately  and  many  cells  are  remained  near  the 
vortex  and  shock  regions  hi  the  other  smooth  regions,  the 
changes  of  the  flow  properties  are  negligible  and  the 
remaining  cells  are  sparsely  distributed.  Fluxes  are 
calculated  only  at  the  remaining  points,  which  enhance  the 
computational  efficiency. 

Figure  7  shows  the  density  contours  of  a  conventional  2nd 
order  accurate  solver,  3rd  order  accurate  solver  and  the 
adaptive  wavelet  scheme  after  the  non-dimensional  time  of 
18.  In  these  figures,  it  is  shown  that  the  contour  of  2nd  order 
accurate  conventional  solver  is  smeared  comparing  with 


that  of  3rd  order  accurate  solver  However,  the  3rd  order 
accurate  adaptive  wavelet  method  can  maintain  the  original 
features  of  a  solution;  the  contours  of  slip  fine,  reflected 
wave  pattern,  etc.  are  clearly  presented  in  the  results  of  the 
3™  order  accurate  solver  and  the  adaptive  wavelet  method. 
Figure  8  shows  the  detailed  density  distributions  of  2nd 
order  accurate  solver.  3rd  order  accurate  solver,  the 
adaptive  wavelet  method  at  y=-5  and  they  are  compared 
with  the  result  of  dense  gnd  system  with  2nd  order  accurate 
solver.  In  the  result  of  2nd  order  accurate  solver,  the  vortex 
and  reflected  wTave  is  dissipated  compared  with  the  3  rd  cider 
accurate  conventional  solver  and  the  wavelet  method. 
However,  by  using  3rd  order  accurate  method,  the  density 
distributions  are  clearly  represented;  the  results  of  3^  cider 
accurate  method  follows  that  of  dense  gnd  system  more 
accurately  though  small  number  of  gnd  points  is  used. 


Figure  6:  An  adaptive  dataset  with  wavelet  level— 3. 


(a)  Conventional  solver  (2nd  order). 


(b)  Conventional  solver  (3  rd  order). 
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(c)  Adaptive  wavelet  method. 


Figure  7:  Comparison  of  density  contour. 


(a)  Density  distribution  at  y=-5. 


(b)  Detailed  density  distribution  at  A  region 


(c)  Detailed  density  distribution  at  B  legion. 
Figure  3:  Comparison  of  density  distribution  at  y=-5. 


4  Conclusions 

Through  this  research,  the  lugher  order  accurate  adaptive 
wavelet  method  is  proposed  for  enhancing  the  prediction 
capability  of  CFD  with  better  efficiency.  For  this  purpose, 
6th  order  accurate  interpolating  polynomial  is  implemented 
to  the  wavelet  decomposition  process.  Also,  the  threshold 
value  is  modified  m  order  to  conserve  the  3rd  order 
accuracy  of  conventional  CFD  schemes.  Hus  lugher  order 
accurate  adaptive  wavelet  method  is  applied  to  one -di¬ 
mensional  shock-sine  wave  interaction  and  two-dimen¬ 
sional  shock-vortex  interaction  problems.  In  these  test,  the 
adaptive  wavelet  method  can  present  much  better  compu¬ 
tational  efficiency  than  a  2nd  order  accurate  scheme  while 
maintaining  the  higher  order  numerical  accuracy  of  a  con¬ 
ventional  CFD  solver. 
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primary  burdens  which  decrease  the  computational  efficiency, 
especially  m  unsteady  calculations.  It  ts  necessary  to  accurately 
compute  the  flow  phenomena  and  efficiently  alleviate  the 
computational  cost  through  the  advancement  of  the  CFD 
algorithms.  The  adaptive  wavelet  method  can  be  a  good 
solution  by  alleviating  the  computational  cost  of  the  spatial 
discretization  [5-7],  The  crucial  idea  of  the  adaptive  wavelet 
method  is  that  the  crucial  locations  m  the  CFD  dataset  such  as 
shock,  vortex  core,  etc.  are  automatically  and  accurately 
searched  by  the  wavelet  transformation.  Then,  the  dataset  is 
automatically  adapted  to  the  local  features  of  a  solution  [5-7]. 
Flux  values  are  computed  only  at  the  remaining  pomts  in  an 
adaptive  dataset,  winch  enhances  the  computational  efficiency. 

In  this  paper,  out  target  is  to  efficiently  and  accurately 
compute  the  Euler  flows  with  shock  discontinuities  For  tins 
purpose,  the  adaptive  wavelet  method  is  applied  to  a  high  order 
accurate  CFD  solver.  For  preserving  the  accuracy  of  the 
original  solver,  higher  order  accurate  wavelet  transformation  is 
constructed  as  follows;  first,  the  6th  order  of  interpolating  poly¬ 
nomial  is  introduced  to  wavelet  decomposition  process. 
Second,  the  threshold  value  is  modified  m  order  to  maintain 
higher  order  accuracy  of  the  solver.  This  higher  order  accurate 
adaptive  wravelet  method  is  fed  mto  the  high  order  accurate 
CFD  solver.  That  is,  a  spatial  discretization  scheme  including 
the  3rd  order  accurate  interpolations  with  e-MLP  hunting 
function  is  performed  at  the  crucial  positions  m  the  CFD 
dataset.  In  the  other  positions,  simple  interpolating  polynomials 
are  used  to  compute  the  residual  values.  To  assess  the  accuracy 
and  efficiency  of  the  adaptive  w  avelet  method,  it  wras  apphed  to 
a  shock-vortex  interaction  problem.  Through  the  application  of 
the  higher  order  accurate  adaptive  wavelet  method,  the 
accuracy  of  the  conventional  solver  is  conserved  and  the  cost  is 
substantially  reduced. 


s'  =  modified  threshold  value 

y  =  ratio  of  specific  heat  (1.4  for  air) 

p  =  density 

THE  OVERALL  PROCEDURE  OF  THE  ADAPTIVE 
WAVELET  METHOD 

The  twro  dimensional  Euler  equation  is  written  as  follows: 


dO  dE  dF 
—  +  —  +  — 
dt  6x  dy 


(1) 


with  O  = 


p 

pu 

,  E  = 

pi 

pr+p 

,  F  = 

P» 

P*v 

pv2+p 

puv 

-per. 

_( peT  +p)u 

( pet+p)v 

and 


a~  1 ,  2  2\ 
- +  —  (u +  V  ), 

yiy- 1)  2V 


wliere  all  properties  and  governing  equations  are  non- 
dimensionalized.  By  the  generalized  coordinate  trans¬ 
formation  Eq.  (2)  is  derived. 


dO  _  dE 
d  r  ~^d% 


dr) 


(2) 


with  Q=$y,  E=j[^Q+4xE+4yF]  and 

F  =  ^j[n,Q+> k-E  +  V7] 


NOMENCLATURE 

a  =  speed  of  sound 

d.j,  =  wavelet  coefficient  at  xt  k 

d"j  =  wravelet  coefficient  at  (i,j)  cell  of  n  tmie  step 
Zj j  =  (p,p,u,v)  at  (/,/)  cell  of  n  time  step 

J2- j  =  interpolated  value  of  2jj 

p  =  pressure 

R”j  =  residual  values  at  (i,f)  cell  of  n  time  step 
R?j  =  interpolated  value  of  Rj'j 

i! ,  v  =  components  of  the  velocity  vector 
Subscript 

i,  j  =  position  of  a  cell 

Superscript 

n  =  current  tune  step 

Greek 

s  =  threshold  value 


FIGURE  1 .  OVERALL  PROCEDURE  OF  THE  ADAPTIVE 
WAVELET  METHOD 
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The  flow  chart  in  figure  1  shows  the  overall  process  of  the 
implementation  of  the  adaptive  wavelet  method  to  a  two- 
dimensional  high  order  accurate  Euler  solver. 

Wavelet  Transformation 

At  first,  assume  a  two-dimensional  dyadic  gndset  as  figure 
2.  Here,  values  at  O  cells  are  saved  to  the  position  m  the 
coarser  level  gndset.  At  the  other  cells  of  □,  A  and  X 
positions,  the  6th  order  of  interpolating  polynomial  is  used  for 
maintaining  the  3rd  order  accuracy  of  a  conventional  solver  as 
shown  m  Eq.  (3). 

n.&W  =^(30"-4.,  -250^+150^, 

+  150&2J-25#Mj+3ef+6j) 

A;  &U  “(36y-4  -25S&-2  +150$,  ^ 

+  150^2-25^4  +  3^6) 

0V,  =  ^(30',j-4  +  3012,2-4  +  20T-2.J-2  -  270,j-2 
- 27012,2-2  +  20V, 2-2  +  30"-4.2  -  270L2 ,2 
X:  +1740-2  +174012.2 -270L4.2+30L6,2 

+  30’-4.2+2  -  2701.2,2+2  + 1 740V+2  + 1 7401.2  2+2 
- 270V.2+2  +  30V,j>2  +  201-2,2+4  -  270’;-4 
-  27012,2+4  +  2014,2+4  +  3012+6  +  30+2,y+6) 


i-2  i  i-+2  i-+4 


o 

□ 

o 

□ 

o 

□ 

o 

A 

X 

A 

X 

A 

X 

A 

O 

□ 

o 

□ 

o 

□ 

o 

A 

X 

A 

X 

A 

X 

A 

O 

□ 

o 

□ 

o 

□ 

o 

A 

X 

A 

X 

A 

X 

A 

O 

□ 

O 

□ 

O 

□ 

O 

FIGURE  2.  TWO-DIMENSIONAL  DYADIC  GRID 


The  difference  values  are  calculated  as  Eq.  (4). 

■  - :  dj+ij  =  Qi+ij  ~  Q!+ij  > 

A:  </+ 


X :  =  G-ij+i  ~  Q-ij+i  - 

After  wavelet  decomposition,  the  wavelet  coefficients  are 
compared  with  the  threshold  value;  if  the  values  are  laiger  than 
the  threshold  value,  the  points  become  remained  m  an  adaptive 
dataset.  Through  these  processes,  the  dataset  follows  the  local 
feature  of  a  CFD  solution.  However,  the  thresholding  process 
inevitably  results  in  the  additional  O(s)  numencal  errors 
because  of  the  data  loss  below  0{s) .  In  order  to  overcome  this 
defect,  the  thresholding  method  is  modified  by  considering  the 
older  of  truncation  errors  as  Eq.  (5)  [7]. 

s'  =  mm[f, max(Ax3 , CFl)  Ax3)] .  (5) 

If  a  wavelet  coefficient  is  smaller  than  s' ,  the  cell  doesn’t 
belong  to  the  dataset.  In  the  other  case,  the  cell  is  mcluded  m 
the  dataset;  an  adaptive  dataset  is  constructed  according  to  the 
local  features  of  the  dataset  with  3rd  order  of  accuracy.  Through 
these  higher  order  wavelet  decomposition  and  modified 
threshold  value.  3rd  order  of  accuracy  of  a  conventional  CFD 
solver  is  guaranteed. 

Flux  Evaluation 

After  the  wravelet  transformation,  flux  values  are  only 
calculated  at  the  mcluded  cells  m  the  adaptive  dataset  with  3rd 
older  of  accuracy.  In  order  to  construct  a  3rd  order  accurate 
spatial  discretization.  e-MLP  limiting  function  is  used  for  the 
oscillation  removal  of  the  3rd  order  interpolation  scheme  [4]. 
Tlie  primary  features  of  e-MLP  are  as  follows;  m  the  CFD,  the 
majority  of  the  regions  are  so  continuous  that  there  is  no  need 
of  applying  the  hunting  functions.  Rather,  the  dissipation  due 
to  hmiting  function  may  damage  the  solution's  accuracy  in  the 
continuous  region.  Thus,  there  is  no  application  of  limiting 
function  m  continuous  regions.  In  the  discontinuous  regions, 
the  regions  can  be  divided  mto  lmear  and  nonlinear 
discontinuous  regions  by  the  types  of  discontinuities.  Then,  by 
the  switching  of  limiting  function  between  TVD  and  MLP. 
TVD  is  used  m  the  linear  discontinuous  regions  and  MLP  is  m 
the  nonlinear  dis -continuous  regions. 

The  bnef  summary  of  3rd  order  accurate  interpolation  with 
e-MLP  is  shown  m  Eq.  (6-9)  [4]. 

(a)  3rd  order  accuracy 


A- 


1  +  2ru 


Pr 


.  1  +  2rRj+l 


(6a) 


,  <*Vi  -  O, 

where  rLi  =  _  . — -  *■- - LJ- 


rRJ-l  = 


~  <£,4 


(4) 


(b)  Continuous  region:  No  limiting  function 
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0£  =  +  0.5/?£A<I>  , 

®R  =  ^i+1  ~  0.5/?*A<I>  - 

(c)  Lmeai-  discontinuous  region: 

TVD  criterion,  <p{r)  =  max(0,min(2,2r)) 

=  <I>,  + 

0.5  max(0. nnn( 2 A<I>? , 2 AO .  ,  pL A<I>?  £  )) 7 
*R  =  <*,+1  “ 

0. 5  inax(0,  nnn(2  AO^ .  2  AOf+1 ,  fiR  AO  f  ))  * 

(d)  Nonlinear  discontinuous  region: 

MLP  limiter  <j>(r)  =  max(0.  min(or,  or)) 

o£  =0,  + 

0.5max(0,  nnn^AO^ ,  aL  AO._^,^£  AO^))  ’ 
*R  =*1+1“ 

0.5  max(0,  min(^AOj+£ ,  aR  AOJ+i ,  ^AO^» 
where  aL  R  are  MLP  coefficients. 


(7a) 

(7b) 


(8a) 


(8b) 


(9a) 
,  (9b) 


After  the  mteipolation  of  primitive  variables  at  the  cell 
interface,  AUSMPW-  scheme  is  apphed  for  the  spatial  discreti¬ 
zation  [8,  9]. 

Residual  Interpolation 

At  unimportant  cells  m  the  dataset,  residual  values  are 
interpolated  by  using  the  6th  order  of  interpolating  polynomial 
as  Eq.  (10). 


□  : 


Klj  =  256OR^j  ~ 25R!L2j  +150^ 


+  150^-251^+3^) 


A: 


=  256(3R'H  ~ 25^-2  + 15°^ 


(10) 


+  150^-25^  +  3^) 


+  2J-4  +  2^-2, /-2  -  VKj-l 

-  27/&2J-2  +  2^4j-2  +  372, ".4 j  ~  2772,"  2j 
x.  +174^+174^-27^+3^, 

+  372,^4  y+2  -  2772,"  2J+2  +  17472^+2  +17472,^J+2 

-  2772,^2  +  372,^2  +  272,".2j+4  -  2772,"J+4 

-  2772^2,2+4  +  272J4J+4  +  3R,nJ+6  +  372,^) 


By  usmg  this  6th  order  of  mterpolatmg  polynomial,  the 
error  caused  by  mteipolation  becomes  smaller  than  the 
truncation  errors  of  the  numerical  schemes.  After  constructing 
residual  distributions  m  the  whole  computational  domain.  3r* 
order  of  Runge-Kutta  tmie  integration  is  performed  [8.  9]. 

NUMERICAL  RESULTS  AND  DISCUSSION 

In  order  to  assess  the  efficiency  of  the  higher  order 
accurate  adaptive  wavelet  method,  it  is  apphed  to  a  shock- 
vortex  inter-action  problem  [10].  The  domain  is  set  as 
-20  <  x  <  5  and  -20  <7  <  10  with  a  400  x  400  gnd.  The 
initial  velocity,  density  and  pressure  distributions  of  a  vortex 
flow  are  presented  m  Eq.  (11)  [10]. 

Tangential  velocity:  ue  =  Mvrexp[(l  -  r2)/2\ , 

Radial  velocity:  ur=  0 , 

Vorticity:  co(r)  =  Mv(2  -  r2)exp[(l  -r2)/ 2],  (11) 


Pressure :  p(r)  =  —  [1  -  - - M~  exp(l  -  r  1  )]r 

/  2 

Density:  p{r)  =  — [1  - ?-^-M2  exp(l  -  r 2 )]1^r_1) . 

The  Mach  number  of  vortex  is  0.39  and  the  initial  vortex  core 
is  located  at  (-5,-5).  The  normal  shock  with  Mach  number  of 
1.29  is  propagated  to  the  vortex  and  s  is  set  as  10°. 

The  overall  efficiency  improvement  according  to  wavelet 
resolution  level  is  summarized  m  Table  1.  In  this  table,  the 
computational  time  with  1  level  of  wavelet  resolution  is  larger 
than  that  of  2nd  order  accurate  conventional  solver.  It  is  because 
Hie  additional  CPU  time  of  wavelet  transformation  is  larger 
than  the  decrease  of  CPU  time  of  flux  evaluation.  By  increasing 
the  wavelet  resolution  level,  the  adaptive  wavelet  method  with 
3rd  order  accuracy  becomes  more  efficient  than  2nd  order 
accurate  solver  as  well  as  3rd  order  accurate  solver;  the 
computational  tmie  of  wavelet  method  became  about  1.7  tmies 
faster  than  the  3rd  order  of  e-MLP  method  at  maxnnum  when 
the  resolution  level  of  wavelet  is  3. 

TABLE  1 .  COMPARISON  OF  CPU  TIMES  IN  SHOCK- 
VORTEX  INTERACTION  PROBLEM 


Level 

CPU  Tmie 

Computational 

Efficiency 

(CPU  sum  fr-ML?/ 

CPU  bsbs Wiy) 

2nd 

order 

e-MLP 

(3"*) 

Wavelet 

1 

1737.6 

2498.2 

1961.9 

1.273 

2 

1537.1 

1.625 

3 

1477.7 

1.691 
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Figure  3  represents  the  adaptive  dataset  with  wavelet 
decomposition  level=3.  The  dataset  follows  the  flow  features 
accurately  and  many  cells  are  remained  near  the  vortex  and 
shock  regions.  In  the  other  smooth  regions,  the  changes  of  the 
flow  properties  are  negligible  and  the  remaining  cells  are 
sparsely  distributed.  By  the  wavelet  transformation,  the  re¬ 
maining  cells  m  the  CFD  dataset  is  about  58%  of  the  whole 
computational  domain  where  fluxes  are  calculated.  In  the  other 
regions,  residual  mterpolation  is  performed.  Here,  the  re- 
maining  regions  are  divided  into  continuous,  linear  dis¬ 
continuous  and  nonlinear  discontinuous  regions;  53%  of  the 
whole  computational  domain  is  continuous  region,  2%  is  linear 
discontinuous  region,  and  3%  is  nonlinear  discontinuous 
region.  Then  m  the  continuous  region,  there  is  no  application  of 
limiting  function.  In  the  hnear  discontinuous  region  and  non¬ 
linear  discontinuous  region.  TVD  and  MLP  limiting  function 
are  used  for  oscillation  removal,  respectively.  Due  to  the  com¬ 
bination  of  wavelet  and  e-MLP.  the  computational  efficiency  is 
substantially  enhanced. 


Guritaiuitr :  52  J*  % 

Lnear  dnsniiiuly  ;  221  % 

■  Naira*  :  2.93% 

FIGURE  3.  Adaptive  dataset  with  wavelet  resolution  level=3 


Figure  4  shows  the  density  contours  of  a  conventional  2nd 
order  accurate  solver,  3rd  order  accurate  solver  and  the  adaptive 
wavelet  scheme  after  the  non-dmiensional  tmie  of  18.  In  these 
figures,  it  is  shown  that  the  contour  of  2nd  order  accurate 
conventional  solver  is  smeared  comparing  with  that  of  3rd  order 
accurate  solver.  Howrever.  the  3rd  order  accurate  adaptive 
wavelet  method  can  maintain  the  original  features  of  a  solution; 
the  contours  of  slip  fine,  reflected  wave  pattern,  etc.  are  clearly 
presented  m  the  results  of  the  3rd  order  accurate  solver  and  the 
adaptive  wravelet  method.  Figure  5  shows  the  detailed  density 
distributions  of  2nd  order  accurate  solver,  3rd  order  accurate 
solver,  the  adaptive  wravelet  method  at  y=-5  and  they  are 
compared  with  the  result  of  dense  grid  system  with  2nd  order 
accurate  solver.  In  the  result  of  2nd  order  accurate  computation. 


the  vortex  and  reflected  wave  is  dissipated  compared  with  the 
3rd  order  accurate  conventional  solver  and  the  wavelet  method. 
Howrever.  by  using  3rd  order  accurate  method,  the  density 
distributions  follows  that  of  dense  grid  system  more  accurately. 


(a)  2nd  order  accurate  computation 


(c)  3rd  order  accurate  e-MLP  with  adaptive  wavelet  method 

FIGURE  4.  Comparison  of  density  contours  with  wavelet  resolution 
level=3 
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X 


(a)  Density  distribution  at 


x 


(b)  Detailed  density  distribution  at  A  region. 


(c)  Detailed  density  distribution  at  B  region 


FIGURE  5  Comparison  of  density  distribution  at  y^5 


CONCLUSIONS 

Through  this  research,  the  higher  order  accurate  adaptive 
wavelet  method  is  proposed  for  enhancing  the  computational 
efficiency  of  e-MLP  method  For  this  purpose,  higher  order 
accurate  wavelet  transformation  is  proposed  including  the 
wavelet  decomposition  with  64  order  accurate  interpolating 
polynomial  and  modified  threshold  value.  This  higher  order 
accurate  adaptive  wavelet  method  is  applied  to  two-dimen¬ 
sional  shock-voiiex  interaction  problems  In  these  test,  the 
adaptive  wavelet  method  can  present  much  better  efficiency 
and  accuracy  compared  with  2“  order  accurate  computation  by 
decreasmg  the  computational  time  of  3m  order  accurate  e-MLP 
method. 
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I,k 


Q"j,k 

Q*j.k 

fl/M 


rv,  v,  >v 


Z 


=  wavelet  coefficient  a|t  ft  /,  A)  cell  of  n  time  step 
=  f(xik),  exact  function  value  at  xiJjr 
=  (p,  p,  it,  v,  tv)  at  C  j,  k)  cell  of  n  time  step 

=  approximated  value  by  interpolating  polynomial 
=  Set  of  original  values  at  Vj 
=  pressure 

=  residual  values  at  (i,  j,  £}  cell  of  n  time  step 
=  interpolated  value  of  Rf  k  given  by  inteipolating  polynomial 
=  components  of  the  velocity  vector 
=  Dyadic  grid  set  with  wavelet  resolution  level 
=  A*11  cell  on  Vi 
=  Integer 


Subsciipt 

ij\  k  =  position  of  a  cell 


Superscript 

n  =  current  time  step 


Greek 

€  =  threshold  value 

e  =  modified  threshold  value 

y  =  ratio  of  specific  heat  (1.4  for  air) 

p  =  density 


I.  Introduction 

RECENTLY  in  aerospace  fields,  Ingh  fidelity  simulations  have  been  strongly  demanded  for  the  accurate  analysis 
of  complicated  flows  related  with  transonic  supers onic  flow  phenomena.  These  demands  emphasize  the 
accurate  and  robust  flow  computation  with  vast  size  of  grids  in  Computational  Fluid  Dynamics  (CFD).  However,  the 
current  level  of  technologies  in  CFD  is  often  insufficient  in  terms  of  accuracy  and  efficiency  ;  a  large  number  of  gild 
pomts  and  higher  order  accurate  algorithms  result  in  a  substantial  burden  in  computation  Sometimes,  the 
consideration  of  computational  efficiency  and  vice  versa  restricts  the  enhancement  of  numerical  accuracy. 

For  solving  these  problems,  many  numerical  approaches  have  been  applied  in  various  CFD  fields.  Among  these 
approaches,  the  adaptive  wavelet  method  has  been  studied  as  very  useful  CFD  tools  for  solutions  algorithms,  feature 
extraction,  and  data  compression,  etc.  In  the  adaptive  wavelet  approach,  wavelet  decomposition  is  performed  to 
calculate  the  wavelet  coefficients  and  some  coefficients  which  are  smaller  than  thresholding  value,  s  are  cut-off  m 
the  dataset.  In  most  cases,  the  w  avelet  coefficients  are  small  in  the  smooth  region  and  large  m  the  rapidly  changing 
regions  such  as  shock,  vortex  and  boundary  layer,  etc  Therefore,  the  crucial  features  m  the  CFD  solutions  are 
automatically  searched  by  wavelet  transformation  and  an  adaptive  dataset  is  constructed  so  as  to  follow  these  local 
features  with  Of)  errors.  Then  high  cost  calculation  such  as  flux  evaluation  is  only  performed  m  these  important 
regions  and  low  cost  interpolation  polynomial  is  applied  in  the  other  legions.  Because  the  major-  parts  of  CFD 
solutions  show  smooth  flow  patterns,  adaptive  wavelet  method  can  present  crucial  advantages  in  the  computational 
efficiency. 

Due  to  this  significant  advantage,  there  have  been  many  researches  about  adaptive  wavelet  methods;  Harten 
presented  an  adaptive  multi-resolution  scheme  for  computing  the  discontinuous  solutions  of  hyperbolic  PDEs.1 
Holmstrdm  proposed  the  algorithm  that  uses  the  inteipolating  wavelet  transformation  to  organize  an  adaptive 
dataset.2  Sjogreen  also  used  a  multi-resolution  scheme  based  on  the  interpolating  wavelet  transformation  to  solve  the 
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compressible  Euler  equations.3  Kang  et  al.  implanted  modified  tines  ho  Id  value  in  adaptive  wavelet  methods, 
applied  the  method  to  two-dimensional  flows,  and  enhance  the  computational  efficiency  of  the  simulation  with  the 
conservation  of  the  numerical  accuracy  of  the  solutions  f  5  D.  Wirasaet  et  al  simulated  three  dimensional  flow 
within  a  differentially  heated  cavity  using  adaptive  wavelet  algorithm  6 

However,  for  more  practical  applications  such  as  three-dimensional  flow  analysis  for  design  optimization  or 
sensitivity  analysis,  the  adaptive  wavelet  algorithm  should  be  simple  to  be  implanted  m  an  original  solver,  enhance 
the  computational  efficiency,  and  presen  e  the  numerical  accuracy  of  the  solver  For  these  purposes.,  the  adaptive 
wavelet  method  proposed  by  Kang  et  al 4  5  is  extended  to  three  dimensional  analysis.  Because  tins  adaptive  wavelet 
method  utilizes  the  order  interpolation  coefficients  as  wavelet  basis,  it  is  very  simple  to  be  united  with  an  original 
solver.  Also,  the  numerical  accuracy  of  an  original  solver  can  be  preserved  by  the  modified  thresholding  value  In 
order  to  assess  the  method,  it  is  applied  to  a  flow  simulation  around  the  GMERA-M6  wing  in  transonic  regime. 
Consequently,  it  is  confirmed  that  the  proposed  three-dimensional  adaptive  wavelet  method  can  enhance  the 
computational  efficiency  of  an  original  solver  with  preservation  of  the  numerical  accuracy  of  the  solver. 

II.  Preliminary  background  of  the  adaptive  wavelet  method 

In  this  research,  die  adaptive  wavelet  method  based  on  the  interpolating  polynomial  is  used."  At  first,  let  us 
assume  a  dyadic  grid  set  as  shown  in  Eq.  (1). 

Vl={^ER:x^=2-ik.keZ}.lEZ.  CD 


The  key  idea  of  the  adaptive  wavelet  method  is  that  die  smoothness  of  flow7  pattern  can  be  easily  determined  by  the 
magnitude  of  difference  value  between  onguial  function  value  and  approximated  value  by  using  neighboring  cells. 
That  is,  the  original  value  can  be  accurately  approximated  by  the  values  at  neighboring  cells  m  the  smooth  region 
but  not  vice  versa  m  rapidly  changing  region.  Here,  interpolating  polynomial  is  used  at  the  odd  numbered  cells  for 
approximation  of  original  function  values  and  then  the  dyadic  dataset  is  decompos  ed  as  Eq.  (2). 


L/r+Ut  -  /tit 

jW+j  =  fu- 

w  here  is  the  order  of  interpolating  polynomial;  the  4th  order  of  mteipolafing  polynomial  is  shown  m  Eq.  (3). 


f  - L  f  +—  f  +  —  f  f 

~  ^  Js.t- i  1 6  ^ ^  16  1 6  jLk~2 ' 


(3) 


After  the  approximation,  the  difference  can  be  computed  as  Eq.  (4). 


dlk  =  J}+i,2fr+i  “  fi+h,7k+i*  Vk  e  Z .  (4) 

In  smooth  region,  the  order  of  difference  value  is  very  small  and  die  cone  spending  cell  can  be  partially  rejected 
m  the  whole  dataset;  an  adaptive  dataset  is  acquired  according  to  the  local  features  of  solution.  And  the  maximum 
error  is  bounded  within  the  order  of  threshold  value  as  Eq.  (5). 

KL=\fv-fi*l<0(e)-  (5) 


III.  Implementation  of  Adaptive  Wavelet  Method  ou  the  Three-Dimensional  Euler  Equations 

In  tins  research,  the  three-dimensional  Euler  equations  are  used  as  the  governing  equations  of  the  flowr  problem.  Hie 
three-dimensional  Euler  equations  are  written  as  Eq.  (6). 
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dQ  6E  dF  dG  n 
dt  dx  dy  dz 
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2  i 

a  i  ?  7  7  •, 

et  = - +  —  (u  +  v  4-  w"  ) , 

7 (y-V  2 


(6) 


where  all  properties  and  governing  equations  are  non-dimensionaiized  By  the  generalized  coordinate  transformation, 
Eq.  (6)  is  rewritten  as  Eq.  (7). 


dQ  3E  dF 
&t  drj 


FA- 

8g 


-Ru,k 


0) 


with  E  =  Ut&  +  ZXE  +  gyF  +  £ZG] ,F  =  Ar,tQ  +  r,xE+nyF  +  r/;G]  and 


G  =jKIQ  +  Z;E  +  CyF+CG]_ 


In  order  to  apply  the  wavelet  procedure  to  CFD  algorithms  of  three  dimensional  Enter  equations,  the  three- 
dimensional  wavelet  decomposition  procedure  is  implemented  to  the  algorithms.  For  preservation  of  the  numerical 
accuracy  of  the  solution  the  modified  threshold  value  is  applied  by  considering  the  grid  spacing.  Hie  flow  chart  in 
Fig.  1  shows  the  overall  procedure  of  implementation  of  the  adaptive  wavelet  method  to  a  three-dimensional  Euler 
solver.  Hie  details  of  all  steps  are  as  follows : 


Yes 


<f  Initial  conditions  > 

i  


Wavelet 

Transformation 


Decomposition 

Thresholding 


Included  in  an  '~ 
L__— — — __A_daptive  Datasel?__. 


Flux  Evaluation 


■  Roe  '  s  Flu*  Difference  Splitting 

■  Koren  Limiter 


Residual 

Residual 

Calculation 

Interpolation 

Time  Integration 
(LU-SGS) 


No 


<  ™  > 

Figure  I  Overall  procedure  of  the  three-dimensional  Euler  solver  with  adaptive  wavelet  method. 
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A.  Wavelet  Transformation 

Assume  a  three-dimensional  dyadic  giidset  with  level  1  as  shown  m  Fig.  2. 


Figure  2  Three-dimensional  dynamic  grids et. 


Here,  values  at  G  cells  are  saved  in  the  coarser  level  gilds  et,  and  4th  order  of  interpolating  polynomial  is 
employed  at  the  other  cells  +  *0  ©)  as  Eq.  (8). 


1  l:  Qi-ljrk  -  —  (-Q?-2j,k  +  PQljJc  +  9Qi+2j,k  -  Q?-4,j:k  ) , 

6y+j,jt  =  +  9%*  +  * "  QIm* >  ■ 

O'  Qi”j,k+I  =  ~(~Qi,j,k-2  +  9&U,k  +  9Qlj,k-2  ~  (?U,k+4 


X: 


.  S+Jj+iJt  “  ~ if +  9Q*jj  +  ^  ;>2,fr  _  0+ti+t-t 


■  fi-2j+JJfc  +  +  9Qi+Zj£  ~  Qi+4j-2jk) 


+ .  QHj-iMi  =  ^-Qlj-2y-2  +  ^ 

“  Q?,j-2tk-4  +  _  Qi!j+4,k-l  ) 

f'j .  0+ij.jt+j  =  —  f-S-zjx-j  +  ^8i",fc  +  -°fl”+2j,Jt+2  ~ 

_  8i'-2jy+4  +  ^8^+2  +  ^fi”+2,j,*  -  Q”+4j,k-2  ) 


(8) 


Qi+lj+lMl  -  —f  -Qi-2j-2,k-2  +  +  ^Sf+2j+2Ffr+2  “  Qt+4J+4A+4 

_  f?^-2j-2F*+4  +  9Q?jMi  +  -  £U;+<*-7  . 

“  fij-2j+4,fc-2  +  +  ^8?-:.j,fr-2  _  0i"+4j-2,lt-4 

-  Qi+4j-2,k-2  +  ^0i"+2j;jt  +  ^0^+2jrfr+2  “  8i  2j+4rk+4  ) 
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After  the  interpolation  of  each  cell,  the  difference  values  are  computed  as  Eq.  (9). 


1  1 :  ^f+i j.k  -  £ 

}lljM  -  Q.’llj.k, 

a:  d”J+lk  =  i 

?u+u-eij+u. 

0:  =  £ 

?/,/.*+!  ~  fij.t+1  ’ 

X  1  ^i+lj+lk  = 

fii+lj+l*  “  0+1/+U  * 

(9) 

^ij+U+l  =  0./+U+T  _  Srj+l.Jt-l  - 
1  =  S-l./.jr-l  _  0+Ij.jt+l  > 

©:  ^r+IJ+U+1  =  0+1, j+U-1  “  6i+\f+\k+1  ■ 


Then  the  difference  values  are  compared  with  the  threshold  value,  e:  if  the  values  are  larger  tlian  the  threshold  value, 
the  points  are  remained  m  an  adaptive  dataset.  If  the  difference  values  are  smaller  than  e.  the  cells  don't  belong  to 
the  adaptive  dataset.  Tliat  is,  an  adaptive  dataset  is  constructed  so  as  to  follow  the  local  feature  of  a  CFD  solution. 

However  due  to  approximation  and  data  loss  from  thresholding,  some  Oft)  error  is  added  to  the  solution  and 
deteriorate  the  numerical  accuracy  of  a  anginal  solver.  In  order  to  prevent  the  addition  of  this  Oft)  error,  the 
thresholding  method  is  modified  by  considering  the  order  of  truncation  errors  4"  ’ Hie  details  are  as  follows.  In  this 
paper,  because  a  second  order  of  MUSCL  scheme  is  applied  to  the  flux  evaluation,  the  flux  at  ft,  j,  k)  cell  can  be 
wntten  based  on  the  real  value  as  Eq.  (10). 


E^El^+Oite1),  F/  =F”rea,+0(Ay2),  Ol=G^+0( Az2) 


(10) 


According  to  the  definition  of  R*}_k  in  Eq  O)  and  Eq.  (10),  Rf_k  is  derived  as  Eq.  (11) 

(11) 

Also,  Oft)  error  due  to  thresholding  is  added  to  the  solution  as  Eq.  (12). 


E‘l=E"+  o{e),  Fi  =  Fi  +o(s ) ,  G"  =G"+  0{s) .  (12) 


Consequently,  R-'j  k  can  be  written  as  Eq.  (13). 

R"j,}=Sij.k  +  0(Ax2s).  (13) 

Hie  third  type  of  error  is  due  to  residual  interpolation  at  the  excluded  points  of  a  SPR  dataset.  In  this  case,  the 
error  order  is  about  0(Av4)  because  a  4th  order  of  interpolating  polynomial  is  used  for  the  calculation  of  the 
residual  values.  Therefore,  the  error  order  at  the  excluded  pomts  of  a  SPR  dataset  is  finally  deiived  as  Eq.  (14). 

«Pj  =  Kj-real  +  0(Al4)  +  0{bx2s)  .  (14) 

__  J 

From  Eq.  (14),  0(e)<0(Ax~)  condition  has  to  be  satisfied  to  conserve  the  second  order  spatial  accuracy. 
Therefore,  for  switching  the  threshold  value  associated  to  grid  spacmg.  the  threshold  value  need  to  be  modified  as 
Eq.(15)_ 


e'  =  min  (s,  Ax') 
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In  this  study,  e  is  set  to  1*10  5,  so  P  =  mm  (10'5,A\:2) . 

B.  Flux  Evaluation 

After  tiie  wavelet  transformation^  flux  values  are  only  calculated  at  the  included  cells  in  the  adaptive  dataset.  In 
this  research,  Finite  Volume  Method  (FVM)  is  used  to  calculate  the  numerical  fluxes;  flux  values  are  then  computed 
on  all  sides  of  the  cell  with  2nd  order  accuracy 

C.  Residual  Interpolation  &  Time  Integration 

At  unimportant  cells  which  doesn't  belong  to  the  adaptive  dataset,  residual  values  are  interpolated  by  4th  order  of 
interpolating  polynomial  as  Eq.  (16). 


□:  =~{  St-2j^  +  9Kj.k  +  ~  XUlkh 

A:  X&u  =  +  9Ki.t  +  9Ki+V  -  Knik), 

O:  R*M1  =  ~ (  8*0-2  +  +  9B”jM2  -  R''JM4  ) , 

.  .  Rihj+l.k  =  —  (~Ri>-2j-2]k  +  91 R”]k  +  9R?+2j+2ji  ~  P?+4,]+4,k 

X .  3z  f 

~  Ri-2j+4,k  +  9Rlj+2£  +  -  ^i+4j-2,k) 

+  .  R-iJ-IMl  =  (~^U~2,k-2  +  9^ijk  +  QR’iJ+Zk+l  ~  ^lj-4rk+4 

-  K HM4  +  9Km2  +  -  KiHk-2  > 

Q .  =  j ~ 

“  K-2JM4  +  9RUM2  +  9Ri+2j,k  ~  Rj+4j]k-2) 

& illj+lk+l  =  ~^i-2j-2,k-2  +  9R*j,k  +  9Ri~2j+2M2  ~  ^i+4j+4.k+4 

Q-  “  K-2j-2,k+4  +  +  9^i+2J+2rk  ~  K+4j+4£-2 

~  K-2J+U-2  +  9^iJ+2,k  +  9K+2,jM2  ~  K-4rJ-2M4 
~  ^i+4j-2,k-2  +  9^i+2j,k  +  9Rtj+2M2  ~  K-2J+4M4  ) 
After  residual  interpolation.  time  integration  such  as  LU-SGS  method  is  applied 

IV.  Numerical  Tests  and  Discussion 


(16) 


For  the  assessment  of  the  computational  efficiency  and  numerical  accuracy  of  the  three-dimensional  adaptive 
wavelet  algorithm,  it  is  applied  to  the  computation  of  flow  field  around  the  ONERA-M6  wing  with  Mach  number  of 
0.8395  and  angle  of  attack  of  3  .06°/  After  the  computations,  the  adaptive  dataset  of  all  level  of  resolution  are  shown 
at  Fig.  3.  These  figures  show  that  the  dataset  follows  the  local  features  of  the  solution;  many  cells  remain  near  the 
wing  and  downstream,  hi  other  smooth  regions,  the  remaining  cells  are  sparsely  distributed  because  the  changes  m 
the  flow  properties  are  negligible. 
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(a)  Level  1 


(b)  Level  2 


(c)  Level  3 

Figure  3.  Adaptive  datasets  of  OXERA-M6  wing  flow  problem. 

The  overall  efficiency  improvement  and  the  L-2  error  according  to  the  wavelet  resolution  level  are  summarized  in 
Table  1  Here,  the  computation  speed  of  adaptive  wavelet  method  becomes  about  3  times  faster  than  the  reference 
CFD  method.  And  the  L-2  errors  betw  een  The  results  of  reference  method  and  adaptive  wavelet  methods  are  only 
about  0(10'  ).  Throughout  tills  implementation,  it  is  substantiated  that  the  adaptive  wavelet  method  can  enhance  the 
computational  efficiency  of  the  steady  state  three  dimensional  Euler  flowr  computation  while  numerical  accuracy  of 
the  reference  CFD  solver  is  preserved. 


Table  1.  Results  of  efficiency  improvements  and  L2  error 


Iteration 

SPRPt. 

CPU  Time 

Tune  Ratio 

L:  Error 

Reference  Solver 

10764 

■ 

29329.3 

- 

■ 

Adaptive 

Wavelet 

1  Level 

3206 

64.61% 

9192.5 

3  19 

2.51  x  10'7 

2  Level 

3267 

64.57% 

9438.4 

3.11 

3.14  x  10'7 

3  Level 

3205 

64.74% 

9251.8 

3.17 

3  .49  x  10  '7 

Fig.  4  shows  the  pressure  contours  on  the  upper  surfaces  of  ONERA-M6  wings  computed  by  reference  CFD 
method  and  adaptive  wavelet  method.  In  both  contours,  the  patterns  of  lambda  shock  which  is  the  crucial  feature  of 
ONERA-M6  wing  are  very  similar.  For  more  detailed  comparison  of  the  results,  pressure  distributions  of  reference 
method  and  adaptive  wavelet  method  at  the  positions  of  50%  and  80%  of  wing  span  are  presented  at  Fig.  5  In  the 
figures,  it  is  shown  that  the  positions  of  shock,  shock  strength  and  pressure  distributions,  etc.  are  almost  same. 
Through  the  application,  it  is  verified  that  the  three-dimensional  adaptive  wavelet  method  can  present  much  better 
efficiency  with  the  maintenance  of  numerical  accuracy  of  a  solution. 
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Pressure 
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0.85 
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0.75 

0.7D 

0.65 

0.60 

0.55 

0.50 

0.45 

0.40 


Pressure 


0.55 

0.50 

0.45 


(a)  Reference  solver  (b)  Adaptive  wavelet  -  2  level  of  resolution 

Figure  4.  Pressure  contour  on  the  upper  surface  of  wing. 


0.46 


3 


0.5 


(a)  Pressure  distribution  at  0.5  of  span  (b)  Detailed  pressure  distribution  at  A  region 


(a)  Pressure  distribution  at  0.8  of  span  (b)  Detailed  pressure  distribution  at  B  region 

Figure  5.  Pressure  distribution  on  the  surface  of  wing. 
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V.  Conclusion 

In  this  research,  the  three- dimensional  adaptive  wavelet  method  is  presented  for  enhancing  the  computational 
efficiency  of  three-dimensional  Euler  flow  simulations  in  transonic  regime  To  this  end,  the  previous  adaptive 
algorithm  proposed  by  Kang  et  al.  is  extended  to  three-dimensional  adaptive  wavelet  properly:  three-dimensional 
wavelet  decomposition  process  and  modified  threshold  value  based  on  the  switching  between  s  and  A*  is 
introduced  to  the  method.  The  developed  three-dimensional  adaptive  wavelet  method  is  applied  to  ONERA-M6 
wing  m  transonic  regime.  Through  the  application,  it  is  shown  that  adaptive  wavelet  dataset  are  constructed  so  as  to 
foliowr  local  features  of  solution  automatically.  Also,  it  enhances  the  computational  efficiency  about  3  times  faster 
than  that  of  reference  solver  with  preservation  of  numerical  accuracy  of  a  solution. 
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Abstract:  Based  on  the  multi-resolution  analysis,  efficient  computational 
algorithm  for  three  dimensional  Euler  equations  is  developed.  A  proper 
thresholding  technique  is  applied  to  preserve  the  order  of  the  accuracy  of  the 
original  solver,  hi  order  to  assess  the  efficiency  of  the  proposed  algorithm,  the 
method  is  applied  to  the  computation  of  flow  field  around  ONERA-M6  wing  in 
transonic  regime.  The  computation  time  is  significantly  reduced  by  the  proposed 
multi-resolution  numerical  method,  compared  to  the  original  solver. 

Keywords:  Multi-resolution  analysis.  Wavelets.  Computational  Efficiency.  Three 
Dimensional  Euler  Equations. 

1  Introduction 

New  design  of  various  aeronautical  and  aerospace  transportation  systems  calls  for  accurate 
computations  of  transonic/supersonic  flow,  that  may  encompass  more  physically  conect  and  detailed 
flow  structures.  However,  current  level  of  Computational  Fluid  Dynamics  (CFD)  modeling  and 
simulation  technology  still  falls  short  of  practical  implementation  for  a  full  scale  or  integrated  design. 
For  enhancement  of  prediction  capability  of  CFD  tools,  it  is  indespensible  to  preform  high  fidelity 
simulations  with  concentrated  grid  point  clustering  and  high  order  accurate  algorithms.  However, 
three  dimensional  flow  analysis  and  unsteady  flow  computation  have  to  face  substantial  increase  of 
computational  cost.  Thereby,  it  is  important  to  develop  an  advanced  numerical  method  that  should 
produce  accurate  simulation  output  and  at  the  same  time,  that  alleviate  overall  computational  cost. 

The  computational  efficiency  issue  has  been  reserached  by  numerical  algorithms  via  multi¬ 
resolution  analysis  using  wavelet  basis.  Especially,  the  wavelets  have  been  studied  as  useful  CFD 
tools  for  solution  algorithms.  feature  extraction,  and  data  compression,  etc.  Harten  presented  an 
adaptive  multi-resolution  scheme  for  computing  the  discontinuous  solutions  of  hyperbolic  PDEs  [1]. 
Holmstrom  proposed  the  algorithm  that  uses  the  interpolating  wavelet  transformation  to  organize  an 
adaptive  dataset  [2].  Sjogreen  also  used  a  multi-resolution  scheme  based  on  the  interpolating  wavelet 
transformation  to  solve  the  compressible  Euler  equations  [3].  Lee  extended  ID  supercompact 
wavelets  concept  to  multidimensions  and  demonstrated  big  compression  ratio  for  massive  flow  field 
simulation  data  [4].  Kang  &  Lee  proposed  modified  threshold  value  in  adaptive  wavelet  methods  and 
enhance  the  computational  efficiency  even  with  the  preservation  of  the  numerical  accuracy  of  original 
algorithm  [5.  6.  7]. 

In  the  multi-resolution  analysis  using  wavelets,  wavelet  decomposition  is  performed  to  compute 
the  scale  and  wavelet  coefficients.  In  most  cases,  quite  many  wavelet  coefficients  are  smaller  than  a 
certain  value.  Even  with  cutting-off  these  small  wavelets  coefficients,  the  reconstructed  data  is  not 
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much  different  from  the  original  dataset.  By  using  efficient  wavelet  data  representation,  a  variety  of 
solution  algorithms  have  been  developed  based  on  multi-resolution  analysis.  However,  practical 
applications  such  as  three-dimensional  wavelet  analysis  for  Euler  and  Navier-Stokes  equations  are 
pretty  rare,  compared  to  the  efforts  on  fundamental  wavelet  analysis.  D.  Wirasaet  et  al.  extended  the 
adaptive  wavelet  method  to  solve  three-dimensional  flows  [8].  Then,  they  simulated  the  unsteady 
three  dimensional  flow  within  a  differentially  fieated  cavity  using  adaptive  wavelet  algorithm  to 
demonstrate  the  versatility  and  efficiency  of  the  method.  However,  for  more  successful  practical 
implementations,  the  multi-resolution  algorithm  should  be  simple  to  implement  in  original  solver  and 
should  be  robust  in  applying  to  non-uniform  and  distorted  grid  system  with  various  shapes. 

In  this  smdy  our  previous  two-dimensional  adaptive  wavelet  algorithm  [5.  6.  7]  is  extended  to 
three-dimensional  flow  analysis.  To  this  end.  three-dimensional  wavelet  decomposition  process  and 
an  appropriate  modified  thresholding  function  are  proposed.  The  developed  multi-resolution  method 
is  applied  to  flow  field  simulation  around  ONERA-M6  wing  in  transonic  regime.  From  the 
application,  it  is  confirmed  that  the  proposed  algorithm  is  about  3  times  faster  than  original  flow 
solver  and  preserves  order  of  accuracy  of  the  reference  solver.  In  addition,  the  limitation  of  current 
three-dimensional  algorithm  and  on-going  approaches  are  discussed. 


2  Multi-Resolution  Analysis  Using  Wavelets  for  Three-Dimensional 
Euler  Equations 


hi  this  study,  the  three-dimensional  Euler  equations  are  used  as  the  governing  equations  to  solve  the 
flow  problem.  The  three-dimensional  Euler  equations  are  as  follows. 

dQ  dE  dF  dG  „ 

dt  dx  dy  dz 
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where  all  properties  and  governing  equations  are  non-dimensionalized.  By  the  generalized  coordinate 
transformation,  above  equations  are  rewitten  as 

be 


dr  d^  dr/  dg  ; 


with  £>=£,  +  F=±[n,Q  +  nxE  +  nyF  +  r,fi]  and 


G  -  —  [ZtQ  +  CXE  +  CyP  +  CZG]  • 

hi  order  to  apply  the  multi-resolution  analysis  to  CFD  algorithms  of  three  dimensional  Euler 
equations,  the  three-dimensional  wavelet  decomposition  procedure  is  implemented  to  the  algorithms. 
For  preservation  of  the  numerical  accuracy  of  the  solution,  the  modified  threshold  value  is  derived  by 
considering  the  grid  spacing.  Also,  three-dimensional  restriction  method  [5]  is  applied  in  time 
integration  for  enhancement  of  the  convergence  of  the  solution.  The  procedure  of  multi-resolution 
analysis  algorithm  is  as  shown  in  Figure  1 . 
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Initial  conditions 


Figure  1 :  Overall  procedure  of  multi-resolution  analysis  algorithm. 

•  Step  1 :  The  multi-resolution  analysis  is  performed:  by  the  decomposition  and  thresholding, 
computational  domain  is  adapted  to  follow  the  local  feaUires  of  the  solution  such  as  shock, 
vortex,  etc.  Then,  the  positions  of  crucial  features  are  included  in  the  adaptive  dataset. 

•  Step  2:  The  flux  evaluation  is  performed  in  the  including  cells  in  the  adaptive  dataset  with 
conventional  methods. 


•  Step  3:  hi  the  excluded  cells  in  the  adaptive  dataset,  the  residual  is  computed  from  the  results  of 
Step  2  by  residual  interpolation. 

•  Step  4:  The  time  integration  is  carried  out.  Then,  restriction  method  is  applied  to  the  solution: 
negligible  flow  variations  are  restricted  and  the  convergence  of  the  solution  is  enhanced. 


2.1  Wavelet  transformation 

2.1.1  Sparse  Representation  (SPR) 

hi  this  study,  the  multi-resolution  analysis  based  on  the  interpolating  polynomial  is  used  [2].  Assume 
a  dyadic  grid  set  as  shown  in  below  equation. 

Vl={xlkeR:xlk=2-lk.keZ}.leZ  . 

The  dyadic  dataset  is  decomposed  into  two  types  of  coefficient  as  follows. 

1/i+U*  =  fl.k 

fi+i2k+i  =  p/+u*+ 1 (ftJ[  Pa '  fi.k+i •  •  ■  •  •  flk+L ) 

where  P,;  is  the  order  of  interpolating  polynomial  in  an  even  number.  Hie  4th  order  of  interpolating 
polynomial  is 


fl+l,2k+l  -  +  —fl.k  +  —fl,k+ 1  —  ~7zfl,k+2  . 

lo  lo  lo  lo 

Then,  the  wavelet  coefficients  are  defined  as 

^ Ik  =  //+1,2*+1  ~fl+l,2k+ 1’  ^  E  Z  • 

In  the  thresholding  process,  the  wavelet  coefficients  corresponding  to  regions  of  less  importance 
are  partially  rejected,  accordingly  sparse  wTavelet  representation  is  obtained.  Through  the  sparse 
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wavelet  representation,  a  compressed  dataset  is  constructed  with  remaining  data  points.  This 
compressed  dataset  is  the  set  of  function  values  that  survive  after  thresholding  as 


2.1.2  Three  dimensional  wavelet  trasformation 

Assume  a  three-dimensional  dyadic  gridset  with  level  1  as  shown  in  Figure  2. 


oo 


m 

w 


+ 


Figure  2:  Three-dimensional  dyadic  gridset 


Here,  values  at  O  cells  are  saved  in  the  coarser  level  gridset.  and  4th  order  of  interpolating 
polynomial  is  employed  at  the  other  cells  (□,  A,  O,  x,  +  ,  O  ,  ©)  for  maintaining  the  2nd  order 
accuracy  of  a  conventional  solver.  The  4th  order  interpolation  is  as 

=  Jgf-iS-ijt  +  9$.JM  +  . 

Qu+l*  =  ~(-Qu  2J,  +  +  9$J+2J,  ~  Quuj,)* 

O  2  +  9$.i*  +  9&,JM2  -  gjj.u), 

x.  Ql+iJ+lJt  ~~ ~(~Qi-2J-2M  +  9Q?Jlk  +  -90+J,j+2>  ~  Oi+4j+4J, 

~  Q?-2,j+4,k  +  +  9(£+2j,k  ~  @i+4,j-2Jk  ) 

+  .  QHj+lJc+1  =  ~^(~8i*,j-2Jc-2  +  9(£,j,k  +  9@i!j+2Jc+2  ~  &i!j+4Jc+4 
~  @i,j-2,k+4  +  98Fj,k+2  +  9(X,j+2Jc  ~  Q"j+4Jc-2  ) 


Q  .  6i+lJ.k+l  =  "J^(-Qi-2j,k-2  +  9(£,j,k  +  9QF+2,j,k+2  -  Ot+4j,k+4 
~  Q?-2j,k+4  +  9@i,j,k+2  +  9@i+2,j,k  ~  @ii+4j,k-2) 
Qi+lJ+l,k+l  =  —  ( ~Qi-2,j-2Jk-2  +  9Q*,j,k  +  9Qj+2J+2M+2  ~  Q?+4j+4Jf+4 

®:  _  8i'-2j-2,k+4  +  9&!j,k+2  +  9Q?+2j+2,k  ~  8?+4j+4Jk-2 

~  Qi-2j+4Jc-2  +  9&J+2M  +  9Q?+2jM2  ~  Q?+4J-2M4 
~  @i+4j-2,k-2  +  9&2j,k  +  98i,j+2,k+2  ~  ) 
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After  the  interpolation  of  each  cell,  die  wavelet  coefficients  are  computed  as 

di+ij  k  =  Qi+ij.k  ~  Qt+xjjk » 
djj+ik  =  Qij+ik  -  Qi.j+ik  ’ 

^rj.fr+1  =  Qi.j.k+ 1  ~  Qi!j.k+ 1> 

X  :  .j+lk  =  Qi+Xj+lk  ~  Qi+lj+Xk  ’ 

+  :  ^"/+U+ 1  =  0FJ+lfr+1  _  Qj+U+1  * 

^N-1.;\*+1  =  0f+l;\*+1  ”Q+t/.*+l» 

©:  ^f+lJ+U+1  =  S-l ;+U+1  “  0+1/+U+1  ‘ 

Then  the  wavelet  coefficients  are  compared  with  the  threshold  value;  if  the  wavelet  coefficients 
are  larger  than  the  threshold  value,  the  points  are  remained  in  an  adaptive  dataset  so  as  to  follow  the 
local  feature  of  a  CFD  solution.  If  a  wavelet  coefficient  is  smaller  than  e'.  die  cell  doesn't  belong  to 
die  adaptive  dataset.  Here,  in  order  to  prevent  additional  0(e)  due  to  the  data  loss  below  0{e).  the 
thresholding  method  is  modified  by  considering  the  order  of  truncation  errors  as 

£  =  tilings,  Ax2) . 

Therefore  the  adaptive  dataset  is  constructed  with  2nd  order  of  accuracy  and  order  of  accuracy  of  a 
conventional  CFD  solver  is  guaranteed,  hi  diis  study.  £  is  set  to  lxlO'5.  so  £’=min(lxiO'5.  Ax2).  More 
detail  about  modified  thresholding  value  is  hi  reference  4  and  5. 

2.2  Flux  Evaluation 

After  the  wavelet  transformation,  flux  values  are  only  calculated  at  the  included  cells  hi  the 
adaptive  dataset.  Tins  research  uses  the  Finite  Volume  Method  (FVM)  to  calculate  the  numerical 
fluxes;  flux  values  are  then  computed  on  all  sides  of  the  cell  with  2nd  order  accuracy.  Using  these  cell 
interface  values,  spatial  discretization  is  performed  with  Roe's  FDS  method. 


2.3  Residual  Interpolation  &  Time  Integration 

At  unimportant  cells  which  doesn’t  belong  to  the  adaptive  dataset,  residual  values  are  interpolated 
by  4th  order  of  interpolating  polynomial  as  shown  hi  below  equations. 

D:  Ku.k  +  pKhj*  -  4"+W> 

A:  =  T6(~R'^  +  9K’k  +  9K^:k  ~ 

O:  Vtk+t  =  -  Kj*+4)’ 


X: 


.  -Rillj+ljc  ~  +  9Kj.k  +  9K+2J+2M  ~  Ri+4J+4M 


K-2J+4X  +  9Kj+2Jc  +  9^i+2j.k  ^1+4j-2Jt  ) 


_|_  :  ^tj+ljc+l  ~  —(-Xij-2*-2  +  9*Zj.k  +  9^i.j+2jk+2  ~  &ij+4Jr+4 

~  K-2M4  +  9Kj.k+2  +  9Kj+22k  ~  Kj+4M-2) 


Q  .  K+lj.k+l  -  ~(-Ri-2,j,k-2  +  9Kj,k  +  9K+2,j,k+2  -  Ri+4j,k+4 
~  Ri-2,j,k+4  +  9Kj,k+2  +  9^2,j,k  ~  K+4,j,k-2  ) 
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Ri+lJ+lJk+l  ~  ^  (~^i-2,j.2Jc-2  +  9tfj,k  +  9Ri "+2j+2Jc+2  ~  ^iUj+4Jc+4 
“  Pi-Jj-2,k+4  +  9Kj,k+2  +  9K+2,j+2Jc  ~  Pi+4,j+4,k-2 
~  & i-2j+4Jc-2  +  9Kj+2Jc  +  9K+2j,k+2  ~  ^i+4,j-2Jk+4 
~  K+4,j-2M-2  +  9tf+2j,k  +  9tfj+2Jf+2  ~  Ri-2,j+4Jc+4  ) 

Through  the  4th  order  of  iiiteipolatiiig  polynomial,  the  error  caused  by  interpolation  becomes 
smaller  than  the  truncation  errors  of  the  numerical  schemes.  After  residual  interpolation,  time 
integration  (LU-SGS)  is  performed  in  the  whole  computation  domain. 

3  Numerical  Tests  and  Discussion 

For  the  assessment  of  the  computational  efficiency  and  numerical  accuracy  of  the  three- 
dimensional  multi-resolution  algorithm,  it  is  applied  to  the  simulation  of  flow  field  around  the 
ONERA-M6  wing  with  Mach  number  of  0.8395  and  angle  of  attack  of  3.06°  [9].  After  the 
computations,  the  adaptive  dataset  of  all  level  resolutions  are  shown  in  Figure  3.  This  figure  shows 
that  the  dataset  follows  the  local  features  of  the  solution;  many  cells  remain  near  the  wing  and 
downstream,  hi  other  smooth  regions,  the  remaining  cells  are  sparsely  distributed,  because  the 
changes  in  the  flow  properties  are  negligible. 


(a)  Level  1  (b)  Level  2  (c)  Level  3 

Figure  3:  Adaptive  dataset  of  three-dimensional  gridset. 

The  overall  efficiency  improvement  and  the  Li  error  according  to  the  wavelet  resolution  level  are 
summarized  in  Table  1.  Here,  the  computation  speed  of  adaptive  wavelet  method  becomes  about  3 
times  faster  than  the  reference  CFD  method,  and  in  this  problem  1  level  resolution  of  wavelet  is  die 
fastest.  The  L:  errors  between  the  results  of  reference  method  and  adaptive  wavelet  methods  are  only 
about  0(1  O'7). 


Table  1 .  Results  of  efficiency  improvements  and  L:  error 


Iteration 

SPRPt. 

CPU  Tune 

Tune  Ratio 

L2  Error 

Reference  Solver 

10764 

• 

29329.3 

• 

• 

Multi- 

Resolution 

Analysts 

1  Level 

3206 

64.61% 

9192.5 

3.19 

2.51  x  10’7 

2  Level 

3267 

64.57% 

9438.4 

3.11 

3.14  x  10*7 

3  Level 

3205 

64.74% 

9351.8 

3.14 

3  .49  x  10'7 
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(a)  Reference  solver 


(b)  Multi-resolution  analysis  -  2  level  of  resolution 


Figure  4:  Pressure  contour  on  the  upper  surface  of  the  wing. 


(a)  Pressure  distribution  at  0.5  of  span 


(b)  Detailed  pressure  distribution  at  A  region 


(a)  Pressure  distribution  at  0.8  of  span  (b)  Detailed  pressure  distribution  at  B  region 

Figure  5:  Pressure  distribution  on  the  surface  of  the  wing. 


Figure  4  shows  the  pressure  contours  on  the  upper  surfaces  of  ONERA-M6  wings  computed  by 
die  reference  CFD  solver  and  die  multi-resolution  analysis.  In  both  contours,  the  patterns  of  lambda 
shock  which  is  the  crucial  feature  of  ONERA-M6  wing  are  very  similar.  For  more  detailed 
comparison  of  the  results,  pressure  distributions  of  the  two  solvers  at  the  positions  of  50%  and  80%  of 
wing  span  are  presented  at  Figure  5.  hi  die  figures,  it  is  shown  diat  the  positions  of  shock,  shock 
strength  and  pressure  distributions  are  almost  same.  Through  the  application,  it  is  substantiated  that 


7 


7 


Published  Paper  in  ICCFD7 


the  three-dimensional  multi-resolution  analysis  can  present  much  better  efficiency  while  numerical 
accuracy  of  the  reference  CFD  solver  is  preserv  ed. 


(a)  Downstream  of  wmg 


(b)  Near  the  wmg 


(e)  Upstream  of  wmg 


Figure  6:  Adaptive  dataset  of  2  level  resolution  at  specific  flow  field. 


Figure  7:  Adaptive  dataset  at  specific  flow  field. 


However  many  data  points  in  the  upstream  and  far  regions  from  the  wing  are  still  remained  as 
shown  in  Figure  6  and  7.  Because  the  interpolating  coefficients  of  this  study  are  derived  from  the 
uniform  gridset.  the  approximation  of  original  value  may  be  insurfficient  where  grid  points  are 
clustered  or  distorted.  In  order  to  solve  this  problem  new  concepts  of  interpolation  or  thresholding 
strategies  are  necessary  such  as  the  reconstruction  of  the  supporting  points  for  interpolation,  the 
combination  of  various  order  of  polynomials,  regional  variable  threshold  values,  and  etc.  based  on 
cell  aspect  ratio  and  skewness.  Howrever.  they  should  originate  more  computational  cost  and  memory. 
Thus  our  on-going  researches  are  the  combination  of  various  order  control  of  polynomial  and  regional 
variable  threshold  values  for  the  more  enhancement  of  computational  efficiency. 

4  Conclusion  and  Future  Work 

hi  this  study,  the  three  dimensional  multi-resolution  algorithm  is  developed  to  enhance  computational 
efficiency  of  three-dimensional  flowr  field  simulation.  The  developed  three  dimensional  multi¬ 
resolution  algorithm  is  applied  to  ONERA-M6  wing  in  transonic  regime.  For  the  application, 
adaptive  wravelet  dataset  are  constructed  so  as  to  capture  local  features  of  solution  automatically. 
However  many  data  points,  including  in  the  upstream  region  of  the  wing,  are  not  discarded  since  the 
current  wavelet  procedure  does  not  reflect  the  gird  skewness  and  aspect  ratio.  Nevertheless,  the 
proposed  multi-resolution  analysis  significantly  enhances  the  computational  efficiency  about  3  times 
faster  than  reference  solver  even  with  preservation  of  the  accuracy  of  numerical  solution.  If  newT 
ideas  such  as  combination  of  interpolation  order  control  method  and  regional  variable  threshoing 
techniequ  are  implemented  based  on  grid  skewness,  it  is  expected  that  computational  efficiency  is 
much  more  enhanced. 
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